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Introduction 



Phase transition is a common term for a wide range of phenomena generally described as 
a transition between different states of matter with the variation of one or more physical 
parameters of the system. The transition is accompanied by an abrupt change of some of 
its physical quantities or its derivatives, whereas the relevant physical magnitudes change 
continuously within a given phase. The simplest example of a phase transition is the 
melting of a crystal to form a fluid when its temperature is increased, which produces a 
discontinuous behavior of its density and some other physical properties. Historically, the 
first classification of the phase transitions was given by P. Ehrenfest [Ehr33] and relies 
on a definition of a phase as a state with the minimum thermodynamic free energy. The 
first-order transition in this framework is a transition with an abrupt change of the first 
derivative of the system's free energy with respect to a certain parameter. The second-order 
transitions are those when the first derivative has a cusp when the parameter is changed, 
that is when a finite discontinuity appears in the second derivative of the free energy. 
Ehrenfest's original proposal was later extended to the cases of infinite discontinuities 
of physical parameters. The higher-order transitions are defined similarly, as the ones 
possessing a discontinuity in n-th (n > 1) derivative of the free energy with respect to 
the parameter. Landau theory [LL80] describes the second-order phase transitions as 
a result of a symmetry breaking, with a rapid change of a so-called order parameter, 
characterizing the symmetry properties of a phase. Well-known examples of second-order 
phase transitions are the transitions between a normal fluid and a superfluid, with the 
superfluid fraction being the order parameter, or the ferromagnetic-paramagnetic transition 
with the magnetization as order parameter. 

Quantum phase transitions are a broad subclass of these phenomena related to quantum 
matter, most generally described as a transition between phases at zero or low enough 
temperature, where quantum effects play an important role. The profound difference of 
quantum phase transitions from the classic phase transitions lay in the absence of entropy 
due to the Nernst heat theorem [Ner07]. A classical description of a zero-temperature 
system can prescribe only one phase (an ideal crystal), whilst a quantum system is capable 
to undergo a transition, but only with the change of a certain non-thermal parameter, 
as for instance its density. The role of entropy in classical systems is played in quantum 
phases by quantum fluctuations. One of the first experimental evidences of a quantum 
phase transition was the solidification to hep solid 4 He at low temperatures with a growth 
of pressure, made by Keesom [Kee42]. The recent advances in methods of manipulation 
of ultracold matter, especially in the topics of cooling and trapping of atoms [Chu98, 
CT98, Phi98] and Feshbach resonances [TTHK99, CGJT10], demonstrated possibilities to 
produce systems with unique and highly tunable interparticle potentials [CTGOll]. The 
tunability of the interaction in terms of non-thermal parameters, which was achieved in 
a number of experiments, plays a key role for quantum phase transitions. One of the 
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first theoretical proposals for a quantum phase transition was the bosonic superfluid-Mott 
insulator transition, based on the Bose-Hubbard model [FWGF89, JBC + 98], that was 
finally experimentally confirmed in the work of Greiner et al. [GME + 02] and a number 
of subsequent experimental set-ups [OTF + 01, TOPK06, TCF+11]. Recently, the phase 
diagram and essential thermodynamics of the three-dimensional Bose-Hubbard model was 
obtained in quantum Monte Carlo simulations by Capogrosso-Sansone et al. [CSPS07]. 

The energy of a quantum system, described by the Schrodinger equation in a state of a 
certain symmetry can be obtained with the help of Quantum Monte Carlo methods. From 
the equations of state, corresponding to different phases, one can find the pressure as a 
function of the relevant parameters. The double-tangent Maxwell construction, based on 
the equality of pressures and chemical potential along the transition line, allows to obtain 
the first-order transition point and the width of the transition zone. 

Quantum Monte Carlo (QMC) techniques are ah initio quantum calculation algorithms 
that might provide deep insight into the design of quantum matter, with a capability 
to describe a multitude of relevant properties and phenomena of the system. Among 
them the possibility to locate quantum and temperature phase transitions, and to quantify 
correlations in the system (e.g. pair correction function, structure factors, and even non- 
local properties, such as superfluid fraction and Bose-Einstein condensate). Bose-Einstein 
condensation (BEC), i.e., a macroscopic occupation of the zero-momentum quantum state 
of a system, despite being proposed by A. Einstein and S. N. Bose in the mid-twenties of 
the previous century [Ein24], [Ein25], used to be considered for many decades more as a 
mathematical abstraction than an achievable state of matter. The superfluid properties of 
4 He at low temperatures, found in the experiments of Kapitza, Allen and Misener [Kap38, 
AM38] are believed to be related to the presence of a Bose-Einstein condensate. The BEC- 
like phase transitions have also been observed in excitonic systems [LW93]. Long-lasting 
efforts of numerous experimental groups to actually observe a signature of a condensed 
phase in ultracold gases finally gave a positive result: in 1995 Bose-Einstein condensate 
was found by E. Cornell and C. Wieman [AEM + 95] in gaseous 87 Rb and later the same 
year in the other alkali vapours of 23 Na (W. Ketterle et al. [DMA+95]) and 7 Li (R. Hulet 
et al. [BSTH95]). The experimental set-ups to produce a condensate are generally quite 
complex, partially due to the strict temperature requirements (T ~ 10~ 7 K). From the 
theoretical sight, QMC simulations can yield accurate predictions about the properties of 
Bose-Einstein condensate, provided the interaction in the system is known. 

Let us explain in more detail the techniques, challenges and results that we present 
in this thesis. We are usually concerned with the properties of a bulk system in its ther- 
modynamic limit, but its QMC description of a bulk is generally performed with periodic 
boundary conditions (p.b.c.) applied to finite size system. Therefore any Quantum Monte 
Carlo method yields results with a certain error, related to the size of a simulated system 
and one has to study the limit N p — > oo, with N p the number of particles. The properties 
of a system in the thermodynamic limit are therefore found out by extrapolating the data 
for limited system sizes to infinity. In the condensed systems that we consider, the conver- 
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gence in the energy goes as 1/N P . The dependence of a certain physical quantity is then 
found for a set of different numbers of particles, and the result is extrapolated to infinity. 
Nevertheless the lowest number of particles, for which the asymptotic 1/N P behavior is 
reached within an acceptable precision is a priori unknown. In practice, a number of probe 
simulations is performed in order to observe the needed linear dependence, but sometimes 
the accessible system sizes are too small to obtain the thermodynamic limit correctly. In 
some cases this problem can be greatly relieved by using the Ewald summation technique. 
In the framework of the Ewald method the potential energy is calculated for the infinite 
system, consisting of periodically replicated copies of the original simulation box. Although 
the method makes the calculation more time-consuming, its use shows that in particular 
systems, especially when the long-range interactions are present, the 1/N P dependence of 
the correction is substantially reduced, making the extrapolation to the thermodynamic 
limit possible. 

The Ewald technique is a well-known simulation tool, often used with some modifica- 
tions in a number of applications, involving molecular dynamics, Monte Carlo and other 
algorithms. Despite its popularity, the scope of its utilization is mostly limited to the 
summation of the Coulomb interactions, where its use is essential due to the long-range 
nature of the potential. Nonetheless, conceptually the Ewald technique may be applied to 
a broad variety of various pairwise potentials, for example, of the generic power-law type 
l/|r| fc . In this Thesis we present a detailed step-by-step derivation of the Ewald sums for 
power-law interaction potentials and for all the terms we give explicit formulas, ready to 
be used straightforwardly in actual simulations. The derived expressions have been used 
in the simulations of systems consisting of Rydberg atoms and particles interacting via the 
Yukawa interaction potential. 

The Yukawa potential has been used in the past as the simplest model interaction in 
atomic nuclei, in dusty plasmas and other systems, but recently this interaction appeared 
in the field of ultracold gases. In recent experiments [TVA + 08, HTY + 11] ultracold systems 
made up of two kinds of fermions, one heavy and another light, have been realized in actual 
setups, with an effective cooling achieved by means of an additional bosonic component. 
Theoretic treatment of quasi-two-dimensional systems with this kind of fermionic mixtures 
has been done in Ref. [PAP + 07]. It is argued that the effective interaction potential between 
light-heavy pairs of fermions is of Yukawa (screened Coulomb) type, with a feasibility 
of reaching a gas-crystal phase transition in two-dimensional geometry. In the present 
Thesis we extend the study of crystallization to a fully three-dimensional case at zero 
temperature. A similar problem was pursued before [CCK76], but unfortunately this 
problem was not solved entirely, and an approximate Lindemann criterion was used instead 
of full scale quantum simulations. We find by means of the diffusion Monte Carlo method 
the exact phase diagram as a function of the relevant parameters, that is density and mass 
ratio between the two Fermi species. The obtained diagram provides valuable information 
on the minimum requirements for the mass ratio to achieve a phase transition in actual 
experiments. Thanks to advances in the field of optical lattices there arises a possibility to 
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produce particle mixtures with extremely high ratios of effective masses. In this Thesis we 
argue that certain existing setups, involving optical lattices, allow to increase the effective 
mass ratio enough for potentially reaching the crystallization. 

In the last decade, there is a new wave of interest in ultracold systems consisting of 
Rydberg atoms, and a number of interesting experiments has been performed [HRB + 07, 
LWK + 09, HLW + 11]. A Rydberg atom is a neutral atom with a single electron excited 
to a high orbital. The important properties of this quantum object are its simplicity 
and similarity to a hydrogen atom. Furthermore, its unique properties of possessing very 
strong and controllable interactions over long distances, together with the novel techniques 
of ultracold atom manipulation, attracted a great deal of attention due to a prospectively 
rich behavior of mixtures of excited and unexcited atoms. In a typical experiment a 
trapped cloud of cold atoms is exposed to a laser field, exciting a small fraction of atoms 
to a particular Rydberg state. These excited states interact in a much stronger way among 
themselves than with the unexcited background. The possibility of tuning, turning on, and 
off, the large magnitude of the forces, as well as a number of other advantageous properties, 
suggest its use as a quantum gate, that is, a basic element of quantum circuits. Currently, 
there is a wide range of proposals for physical systems to realize quantum information 
processing units: trapped ions [BW08], linear optics [KMN + 07], superconductors [CW08], 
quantum dots [LWS + 03], and so on. The one, based on systems made up of Rydberg 
atoms is unique in terms of the range and the amplitude of the interaction, the working 
frequency and other advantageous properties [SWM10]. The basic principles of a trapped 
Rydberg system as a quantum gate stem from the idea of a so-called Rydberg blockade, 
that is when a single excited atom shifts the energy levels of the nearby unexcited atoms 
out of the resonance with the driving laser pulse. Further excitations, injected into the 
cloud, can bring a macroscopic fraction of the cloud to a blockaded mode, allowing for 
a partial or complete saturation. In actual experiments, the fraction of unexcited atoms 
permits over 10 3 excitations before the suppression of the new ones appear [HRB + 07, 
TFS + 04, SRLA + 04]. The actual physical phase of the excited atoms cannot be accessed 
directly in the reported experiments, although their arrangement is considered as a relevant 
information, both as a standalone physical problem and for the implementation of the 
quantum gate. A direct observation of a quantum phase transition and a presence of long- 
range ordering is argued to be a feasible task in similar systems [WLPB08, LWK+09]. In 
the field of quantum computations, Pohl et al. [PDL10] proposed that the presence of a 
crystal-like ordering could provide a better control over the quantum states. An insight to 
the spatial ordering in a cloud of Rydberg atoms may also shed light to the phenomenon 
of the so-called "antiblockade" [APPR07]. 

As mentioned before, the behavior of an ultracold mixture of excited Rydberg atoms 
and unexcited background is profoundly rich and complex. It can also depend substantially 
on particular experimental conditions, like the cloud geometry, laser field properties, etc. 
We perform a study of a model system, in which we neglect the interactions related to 
the unexcited background, and using the pairwise repulsive van der Waals C 6 /|r| 6 for the 
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excited atoms. The general aim of this study is to fully understand the phase diagram of 
the system. A perspective comparison with future experimental results can demonstrate, 
how well the properties of the system can be derived from this simple model. Since the 
number of Rydberg atoms, typically present in the current experimental works, is of the 
order of thousands or greater, in our simulations we look for all the relevant results in the 
thermodynamic limit. There is also a variety of possible crystal packings, which might be 
realized in the solid phase, hence we give a discussion on which of them are energetically 
preferable. 

Another physically relevant system, considered in the Thesis, is bulk molecular para- 
hydrogen (p-H 2 ). This system in the quantum regime (at low temperature) was proposed 
theoretically as a possible candidate for superfluidity, but it crystallizes at the tempera- 
ture substantially higher than transition temperature, making it impossible to observe a 
transition to the superffuid phase. In this work, our Group has studied a metastable non 
crystalline phase of bulk p-H 2 by means of Quantum Monte Carlo methods in order to find 
out the temperature at which this system still contains a noticeable superffuid fraction. 
The ultimate goal that our Group pursued, was to frustrate the formation of the crystal in 
the simulated system and to calculate the temperature dependence of the one-body density 
matrix and of the superfluid fraction. I present the study of the limit of zero temperature 
using the diffusion Monte Carlo method. Results for the energy, condensate fraction, and 
structure of the metastable liquid phase at T = are reported and compared with the ones 
of the stable solid phase. The simulation at zero temperature is used by our Group as a 
starting point for the simulation of the system at low temperatures by using Path Integral 
Monte Carlo technique. 

The structure of the Thesis is as follows. 

In Chapter 1 we discuss the analytical approaches and approximations used in the 
subsequent Chapters; also we describe the general concepts of the two-particle scattering 
problem as a tool to construct Jastrow terms in trial wave functions. Chapter 2 explains in 
details the Quantum Monte Carlo methods employed in our calculations from the theoreti- 
cal and practical points of view. In Chapter 3 we explain the Ewald summation technique, 
applied to a power-law l/|r| fc interaction potential, and a generic approach to obtain the 
Ewald terms. The obtained expressions of this analytic work are implemented into sim- 
ulations of different physically relevant systems (Rydberg atoms and Yukawa particles). 
Chapter 5 is devoted to the modelling of a system, governed by the model potential be- 
tween Rydberg atoms l/|r| 6 . The phase diagram of the system is obtained for a relevant 
range of densities and temperatures, combining quantum simulations at low temperature 
and classical treatment at higher temperature. A special attention is paid to the clas- 
sical description of this system, composed of Rydberg atoms, and its comparison to the 
quantum system. In Chapter 4 we present the simulation of a system with the Yukawa 
interaction potential. The following Chapter 6 presents the results of the Quantum Monte 
Carlo simulations of molecular para-hydrogen at zero and finite temperatures, performed 
in our Group. Conclusions are drawn in Chapter 7. 



Chapter 1 

Tools 



1.1 Introduction 

This Chapter is intended to provide theoretical basis for the following Chapters. The quan- 
tities characterizing properties of a quantum system (correlation functions, static structure 
factor, and so on) are introduced here and are later used in the subsequent chapters. We 
also discuss the two-body scattering problem in a three-dimensional system geometry that 
sheds light to the short-range properties of many-body systems. The two-body scattering 
solution can be used in the development of trial wave functions needed in the Quantum 
Monte Carlo algorithms. 

The structure of the Chapter is the following. 

In Section 1.2 we introduce experimentally relevant magnitudes and functions that are 
present in a quantum system. First of all, we consider the analytic forms of the first 
and second quantization (Sees. 1.2.1, 1.2.2). Special attention is given to the relation 
between correlation functions and mean values of quantum operators. Some correlation 
functions may be greatly simplified in case of a homogeneous system, which is presented 
in Section 1.2.2. The definitions and general comments on static structure factor and 
momentum distribution are drawn in Section 1.2.3. 

Section 1.3 is devoted to the study of two-body scattering processes in three-dimensional 
geometry. The solutions of the two-body scattering problem are provided for a number of 
physically relevant interaction potentials. The main aim of this last Section is to give an 
efficient tool to construct two-body Jastrow terms of the trial wave function for Quantum 
Monte Carlo simulations (for details on QMC methodology see Chapter 2). 



10 



Chapter 1. Tools 



1.2 Correlation functions 
1.2.1 Second quantization form 

A quantum system of identical particles of variable number is generally described with 
the help of annihilation and creation operators. The commonly used notations for the 
auxiliary field operators are \&t(r) for an operator creating a particle in the position r, and 
^(r) for an operator destroying a particle in the same position. By means of the creation 
operator a m for m-th state, that puts a particle to an orbital (p m (r), and the annihilation 
operator, aj^, that removes a particle from the orbital </? m (r), these field operators can be 
easily represented in the following form: 




£</4( r ) S m 



'"m 

m 



If we consider a uniform gaseous system within a volume V, single particle states are evi- 
dently plain waves </? m (r) = e lkmT /\/V. Bosonic operators (1.1) obey commutation relation 
[^(r), ^(r')] = S(r — r'), [^(r), ^/(r')] = 0, while fermionic operators obey commutation 
relations. 

First of all, let us discuss the relation between the correlation functions and the mean 
values of one- and two-body quantum mechanical operators. Let us consider the simplest 
case when the Hamiltonian of the system is a sum of only one-body and two-body operators 

H = FW + F ( - 2) , (1.2) 

where the one-body operator F^ stands for a sum of the one-body terms, and the two- 
body operator F^ is a sum of corresponding two-body terms, depending on r^, Yf. 

N 

= (1.3) 



i=l 

N 



1 iv 

F [2) = lY,fV- (1-4) 

Obvious examples for one-body operators are an external potential field, depending 
only on the particles' coordinates: f^'(r) = V cxt (v), or the kinetic energy: f^(p) = 
p 2 /2m. The first operator is diagonal in coordinate space, while the second one is diagonal 
in momentum representation. A typical example of a two-body operator is a pairwise 
interaction potential, given in coordinate space: /( 2 )(ri,r 2 ) = Vi nt (ri,r 2 ). 

The representation of one- and two-body operators and in terms of field 
operators (see (1.1)) is straightforward: 

= JJ&(T)fV)(T,r')$!(T , )didi f (1.5) 
jM 2 ) = I JJ&(r)&(r')f {2 \r,r')^(r')y(r) drdr' (1.6) 



1.2. Correlation functions 



11 



where the factor 1/2 is introduced to take into account the double summation. 

Up to now, we did not restrict ourselves to local one-body operators (that is those 
satisfying the relation for the quantum averages (rla^^r') = a^(r)5(r — r')), but we also 
consider non-local operators, that is, ones allowed to depend on two arguments = 
a^'(r, r') in the corresponding integral of (1.1). 

Correlation functions can be introduced in terms of the field operators in the following 
way: 

CiO-y) = (^(r^r')) , (1.7) 



C 2 (r,r') = (#t( r )£t( r ')£( r ')£( r )} ? (L8) 

Note that in one-body correlation function (1.7) we consider non-local dependence, and it 
has two arguments. At the same time we consider only local two-body operators, that is 
why we keep two arguments instead of four in (1.6). 

The quantum averages of the operators F^ and may be obtained from p 1 ' and 
p 2 \ when the correlation functions are known: 



(F®) = J fWtfdfar) dr (1.9) 
(F (2) > = IJJ f i2 \r,r')C 2 (r,r')drdr' . (1.10) 

The correlations of the field operator between two distinct points r and r', are char- 
acterized by the one-body correlation function Ci(r, r'). The diagonal component of (1.7) 
r = r' yields the density of the particles p(r) = (W (r)^(r)) = Ci(r, r), hence the sum 
over the diagonal terms, i.e., the trace of the matrix Ci, is equal to the total number 
of particles N = trCi = / Ci(r, r) dr. The two-body correlation function C 2 (r, r') de- 
fines correspondingly the density correlations between the particles at positions r and r', 
respectively. 

It is convenient to introduce dimensionless versions of these functions (1.7) and (1.8): 

t /\ Ci(r, r ) 
Ci(r, r) = ; r= (1.11) 



P(r)\M r ' 

c 2 (r,r') = ^1 (1.12) 
p(r)p(r') 

It can be seen that for bosons the range of any of the functions (1.11-1.12) is [0, 1], and 
the function can be interpreted as a probability to remove a particle from the position r 
and place it to the position r'. The obvious relation ci(r, r) = 1 reflects the fact, that there 
is always a possibility to put the particle back to its initial location. If no Bose-Einstein 
condensate is present, the non-diagonal terms asymptotically vanish in the long range limit 
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ci(r, r') — > 0, |r, r'| — > oo. The function c 2 (r,r') can be understood as a joint probability 
to find one particle in the point r and another one in the point r'. 

Additional information on correlation functions can be found in Refs [Fee67, MahOO, 
Gla63, NG99, GS03]. 



1.2.2 First quantization form 

The physical meaning of the correlation function written in the form of the second quanti- 
zation has been briefly discussed in Section 1.2.1. We will use the Monte Carlo methods in 
order to evaluate averages over the wave function "0(R) of the system. For that one needs 
to represent the averages as integrals of the operators over the wave function ip(R). We 
will look for the mean values of the operators in forms, similar to that of (1.9) and (1.10). 
In the first quantization the expectation value of a one-body operator reads as 

(1) = /^(R)^V(R) dR = gmri,..,!-^ V>(r 1 ,....,r„)dR 
^ ' /|^(R)| 2 dR |^(R)| 2 dR 

~ N /|^(R)| 2 dR ~J a ( r ) c ^ r ) dr > ( L13 ) 

with the notation Ci(r, r') used for the expression 

n( n Ar /^(r, r 2 , ...,r N )ip(r', r 2 , r N ) dr 2 ... dr^y 

Cl(r ' r) = iV lf(r, r.W,,,^^,.^, (L14 > 

An average of a two-body operator (1.10) can be expressed in terms of the two-body 
correlation function (1.8) in the following way: 

N 

M(2 , /f(R^ W (R) dR ^ jV ' (r ' '•v^ 2 ^-- 1 ^^^ ^ dR 

[ ' /|^(R)| 2 dR \ip(R)\ 2 dR * ' } 



N{N- l)ja^{v u v 2 )\^{v u ...,r N )\ 2 dR 1 r (2) 



/|^*(R)| 2 dR 2 



a^(n,r 2 )C 2 (n,r 2 ) dndi&.16) 



with the following expression for the first quantization form of the two-body correlation 
function 

The situation can get easier if we stick to the case of a homogeneous system, as it 
possesses translational symmetry. In the case of one-body operator, the diagonal r = r' 
element of Eq. (1.14) is just a constant, which value is fixed by the density n = N/V, 



1.2. Correlation functions 



13 



Ci(r, r) = n. The non-diagonal elements of the normalized matrix of one-body operator 
in the first quantization (or simply one-body density matrix) read as 

Nfvp*(r,r 2 , ...,r N )vp(0,r 2 , ...,r N ) dr 2 ...dr7v . . 

9i{r) = FT —, 1.18 

n J |^(ri,...,rjv)| 2 dri...drAr 

where n = N/V is the average density of a homogeneous system. The normalized two-body 
density matrix, also called pair distribution function is represented by 

_ N(N - 1) / |^(r, 0, r 3 , -., r^)) 2 dr 3 - dr^y 

92\T) — o TTTt \vTl\ j U- iy J 

n 2 J \ip(Yi , . . . , r jvj | dri...drjv 

The basic properties of the pair distribution function in the zero temperature limit can 
be deduced in the following way. In a gas phase, density- density correlations vanish for large 
interparticle distances, which corresponds to <?2( r ) — > 1 — hence in the thermodynamic 
limit g2{r) asymptotically tends to 1. 

One faces the opposite situation at short distances, where the particle correlations 
are strong, and the value of #2(0) can vary depending on the interaction potential. For 
instance, in case of a repulsive potential #2(0) < 1, on the contrary for a purely attractive 
potential #2(0) > 1, and in the case of a hard-core interaction when the particles are not 
allowed to overlap, #2(0) = 0, when |r| < R core . 

Let us consider the expectation value of a two-body operator (1.10) and see how it can 
be simplified in a homogeneous system: 

Ai2) ) = Y /^(r 1 ,r 2 )a^(|r 1 -r 2 |)dr 1 dr 2 = ^ J g 2 (r) a®(r) dr (1.20) 

In the integration performed above we made use of the mentioned fact, that the operator 
depends only on the distance between particles, allowing to integrate one of the coordinates 
out. 

In the simplest case of a contact delta-potential V int (r) = gS(r) (g stands for a coupling 
constant, defining the interaction strength) the potential energy is simply proportional to 
the value of the pair distribution function at r = 0: 

^f = \gngM (1.21) 
1.2.3 Static structure factor 

By virtue of the field operator (1.1) the momentum distribution is represented as 

n k = (1.22) 
The field operator in momentum space ^ is in fact the Fourier transform of ^(r): 

*(r) = J<?*v m -£J (1 ' 23) 

V Z7T 
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where D stands for a dimensionality of the system. Applying the relation (1.23) to the 
equation (1.22) one finds the following form of the momentum distribution 

nk= L° II eiksCi (I + R; R _ ?) ds dR (L24) 

It can be noticed that in a homogeneous bulk system the center of mass motion can 
be integrated separately, as the dependence on the momentum is defined by the distance 
between the particles, not positions themselves. Performing this integration for a homoge- 
neous system it yields 

n k = n J e~ ifcr #i(r) dr (1.25) 

For example, for a fully Bose-condensed gas, the off-diagonal terms of one-body density 
matrix are constant, gi(r) = n, which results in a singular momentum distribution = 
NS(k). In this case all particles are condensed in k = state. 

Another useful quantity is the dynamic structure factor of the system, which character- 
izes a scattering process, corresponding to the exchange of energy hu: and the momentum 
hk in the scattering event. The dynamic structure factor S(k,uj) can be expressed by 
virtue of the /c-component of the density operator at zero temperature 

TV 

i=i 

in the following form 

i 

where Ui is the frequency of the i-th stationary state. The static structure factor is pro- 
portional to the frequency integral of the dynamic structure factor, that is it characterizes 
the overall probability of scattering of a probe particle with the momentum transfer hk. 
The separate integration over oj gives 

S{k) = NJ S ( k ^)^-^((PkP-k)-\(Pk)\ 2 ) (1-28) 

The latter expression (1.28) can be sampled directly in Quantum Monte Carlo simu- 
lations. The other way to write the static structure factor is to relate it to the two-body 
density matrix by means of the equations (1.8) and (1.26), having in mind the commutation 
properties of the field operator ty(r) 

S(k) = l+U ^ {r2 - Tl)k (C 2 {r 1 , r 2 ) - n{n)n{r 2 )) dr 2 dr 2 (1.29) 

In case of a homogeneous system, the positions of the particles enter the two-body 
density matrix only as an interparticle distance r = |ri — r 2 |, thus the static structure 
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factor can be rewritten in terms of the Fourier transform of the pair distribution (1.12) of 
the system: 

S(k) = 1 + n J e [rk (g 2 (r) - 1) dr (1.30) 

The static structure factor can yield valuable information on the arrangement and the 
order of the system, and its value can be directly accessed in spectroscopy experiments. 



1.3 Scattering problem 
1.3.1 Introduction 

The construction of a trial wave function for a many-body problem is in most cases a very 
complex task since the exact solution is generally unknown (rare exceptions for analytic 
solutions are the ID gas of hard rods and hard points (Tonks-Girardeau gas), where the 
analytic solutions are known). A typical approach to develop a trial wave function therefore 
consists in matching long-range behavior with a two-body solution at short distances. 
In this Section, we will be concerned with the two-body scattering problems in three- 
dimensional systems, whose solutions are then used in many-body calculations presented 
in the following chapters. 



1.3.2 Scattering problem in three-dimensional geometry 
1.3.2.1 General approach 

In this Section, we formulate a generic two-body scattering problem in a 3D geometry. For 
a three-dimensional system the low-density regime of interparticle interaction is supposed 
to be correctly described by two-body collisions. 

Consider two particles with coordinates r x and r 2 and respective masses mi and m 2 
staying close enough to see the process as a two-particle collision. We suppose that the 
system is not confined externally, hence the problem may be treated as translationary 
invariant, with the center of mass moving with a constant speed. Our purpose is to obtain 
the stationary solution p(r!,r 2 ) of the Schrodinger equation 



h 2 



2mi 



-A 



h 2 



ri 



2m 2 



A r2 +y(| ri -r 2 | 



p(r 1 ,r 2 ) = S 12 p(r 1 ,r 2 ) 



(1.31) 



In the center of mass frame, the coordinate variables get separated, thus the representation 
of the Schrodinger equation for the motion of the center of mass 1Z = (miri +m 2 r 2 )/m tot ai 
gets simple 



h 2 

2?Tltotal 



;i.32) 
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with mtotai — m i + m 2 staying for the total mass. The solution of Eq. (1.32) is of a free 
wave form. Omitting the normalization, the solution reads fn(TZ) = exp[ik-jiTZ] with kn, 
being the initial center of mass wave number of the system with the corresponding energy 
£ n = h 2 k^/2m tota i. 

The equation for the position difference r = r x — r 2 contains the pairwise interaction 
potential 

l " A r + V(\T\))f(T) = Sf(T) , (1.33) 



2/i 



with 



p=-Z*p- (1.34) 
mi + m 2 

denoting the reduced mass. When one finds the solutions of Eqs (1.32-1.33), it is possible 
to obtain the needed solution of the scattering problem (1.31) as 

fp(ri,r 3 ) = M(K)f(r) (135) 

&12 = £ + &R. 

We will assume that the energy of the incident particle E is small enough and the 
solution is spherically symmetric /(r) = jf(|r|). Under these prescriptions, and using the 
spherically symmetric representation of the Laplacian A = + 2 jL^ one can conveniently 
rewrite the equation (1.33), introducing the auxiliary function w(r) 

Sir) = ^ (1.36) 
r 

in a way that its analytic form is similar to a one-dimensional Schrodinger equation: 

h 2 

-—w"{r) + V int {r)w{r)=Sw{r) (1.37) 

with the additional requirement of the boundary condition 

w{0) = . (1.38) 

At distances, large compared to the range of the potential, one can neglect VJ n t( r ) term in 
Eq. (1.37) leaving with a free wave differential equation. 

- — w "(r) = Ewir) (1.39) 
The solution for this equation is a plane wave 



w(r) = Asm(k s r + S(k s )), 



(1.40) 
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where 



hk s = yljiE 



(1.41) 



stands for the momentum of the incident particle and 5(k s ) is the scattering phase and A 
is an arbitrary constant. 

Properties of the low-energy scattering depend on a single parameter, namely s-wave 
scattering length a 3 ^- Its value can be obtained from the phase shift S(k s ) as the following 
limit 



If we consider the asymptotic low-momentum limit k — > (slow particles) the scattering 
solution (1.40) may be expanded 



and it is easy to see, that it has a node at a distance a^D- The position of the first node 
of the positive energy scattered solution in the limit of low momenta can be seen as an 
equivalent definition of the scattering length in a three-dimensional system. 

In low-density regime of weekly interacting gas the interparticle distance is large com- 
pared to the range of the potential. Under such conditions the exact shape of the interaction 
potential is not important and the description in terms of the s-wave scattering length is 
universal. 

In the next several sections we will consider the scattering problem for a hard-sphere po- 
tential (1.3.2.2) as a simplest example, and afterwards the same problem for the Yukawa 
potential (1.3.2.3) and the common potential between Rydberg atoms (1.3.2.4) will be 
solved. We will obtain explicitly the expressions for the scattered functions, which are of 
great importance, since in many cases they can give a deep physical insight into properties 
of a many-body system. Indeed, under particular conditions the relation between the cor- 
relation functions and the scattered wave function f(r) can be found. Another important 
point to mention is that often a many-body trial wave function is taken as a product of 
scattered functions f(r). Hence, these calculations are of great importance for their further 
implementation of the Quantum Monte Carlo algorithms. 

1.3.2.2 Scattering on hard sphere potential 

As it was mentioned in Sec. 1.3.2.1, in the limit of low energy collisions the interaction 
potential is characterized exclusively by one parameter, the s-wave scattering length. It 
means that when E — > the scattering on each potential possessing equal scattering 
lengths is the same, and it is said that the scattering becomes universal. If the scattering 



a 3D = - lim 






(1.43) 
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is considered on some repulsive potential, then the easiest choice is the hard sphere (HS) 
potential: 

V^(r) = { +°°' r < ° 3D (1.44) 
[0, r> a 3D 

This potential is controlled by a single parameter, which we denote a 3 o in the definition 
(1.44). It is clear that we can treat this value as a range of the potential. Simultaneously 
it has a meaning of the scattering length, as it was presented in (1.42). It comes out in 
a natural way from the solution of the scattering problem, as by definition the s-wave 
scattering length is the position of the first node of the analytic continuation of a large- 
distance free wave solution. For the hard sphere potential this free-wave solution is valid 
for r > a 3 r> with the position of the node given by |r| = a 3 r>. 

The Schrodinger equation (1.37) can be rewritten as (the reduced mass is \x = ra/2) 

ft 2 

w"(r) + V hs (r)w(r) = £w(r) (1.45) 

m 

A particle is unable to penetrate the potential wall of the hard core thus the solution 
tend to zero for distances below the size of the hard sphere. Notice that the energy in this 
case has only a kinetic component, and the interaction potential does not enter explicitly. 
Nevertheless it affects the boundary condition for the solution. 

[ ™(r) = 0, |r| < a 3D 

\ w"(r) + yCwir) = 0, \r\ > a 3D K J 

with a convenient substitution for the frequency x 2 = ^f. The solution of this differential 
equation (1.46) may be obtained easily. Joining together with (1.36) one gets: 

fir) = I °' |r| < a3D (1 47) 

J{r> \ Asm(k(r - a 3D )) /r, |r| > a 3D ' [i ^ n 

with A for an arbitrary constant. The phase shift is related linearly to the wave vector 
of the incident particle 5(k) = —xa 3 D (1-42), showing explicitly that the range of the 
potential (1.44) has in fact the meaning of the three-dimensional scattering length as we 
said above in the same Section. 



1.3.2.3 Scattering on Yukawa potential 

The solution of the two-body scattering problem, as it was mentioned above, is used 
to construct the trial wave functions for subsequent use in the Quantum Monte Carlo 
simulations. The need of having correctly posed short-distance approximation for the trial 
wave stems from the fact that the majority of physically relevant pairwise interactions 
contain a repulsive core, that can vary in hardness. Since the diffusion Monte Carlo 
method is based on a sampling from the particle distribution ipr^o, & n inaccurate choice 
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of the trial wave function i/jt for small r can result in a substantial growth of the overall 
error of the calculation. This can be manifested by a need of raising the number of walkers 
to reach a convergence as well as by a growth of a common statistical variance. On the 
other hand, the errors in the long-range part are usually easily "corrected" by the DMC 
algorithm, as this range is statistically well represented, and also the actual discrepancy 
of the trial wave and the ground state solution is relatively low. 

The two-body scattering problem for the potential of Yukawa is solved in a similar 
manner, as was used in the solution for a system of hard spheres. The solution of the 
initial Schrodinger equation is considered spherically symmetric and therefore rewritten as 
a one- dimensional equation with a single argument, the interparticle distance r = |ri — r2 1 . 
If we look for the solution in the form convenient for subsequent use as a Bijl-Jastrow term 
in Monte Carlo calculations, then in accordance with Eq. (1.37) 



-— w"(r) + V eXpi - 2r) w(r) = £w(r) . (1.48) 
m r 

The solution of the latter equation for short-range distances can be obtained by approxi- 
mating exp(— 2r) by 1. One possibility is to consider the scattering at a finite energy £ and 
fix the value of £ from continuous matching with the long-range behavior of Bijl-Jastrow 
term. Equation (1.48) can be solved by means of the hypergeo metric functions. A less 
precise description can be obtained by setting £ = 0, still the obtained Bijl-Jastrow term 
is well suited for DMC calculation. The £ = case 

_ ^ w "{r) + 6Xp( ~ 2r) w (r) = (1.49) 
m r 

results in a more simple solution which is a linear combination of the modified Bessel 
functions of the first kind and the second kind with particular square root factors, 

xyr xyr 



with x = yh 2 /(mVo). The first component I\ is finite for r = 0, while the second one K\ 
diverges for r = 0. The arbitrary constant is irrelevant for the QMC algorithms that we 
use, therefore we stick to the following form of the two-body scattering solution 

f(r)= U V ; = x + — r + — r 2 + 0(r 3 ) (1.51) 



A more precise analytic solution can be found by using a higher-order expansion of the 
Yukawa factor exp(— 2r), that is exp(— 2r) « 1 — 2r, however the use of a single first term 
of the expansion is usually enough for practical purposes. 
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The diffusion Monte Carlo study of Yukawa systems, made by Ceperley et al. [CCK78], 
was based on the following form of the Jastrow term of the trial wave function, typically 
used in the nuclear matter calculations 



u[r) 



Ae- Br (l-e- r/D )/r. (1.52) 



It can be shown that the expression exp(— u(r)) can be made coincident in the leading 
terms with the expansion of the solution (1.51) of the two-body scattering problem: 



A(l + 2BD)e~ A / D 
'2D 2 



e -u(r) = e ~A/D + -v- ; -—^ r+ (L53) 



A(3D -AD + Y1ABD - Y2BD 2 + 12AB 2 D 2 - 12B D )e~ ^ 9 , n/ 

H — = r 2 + 0(r s ) 

24D 2 v ' 

One can note that Jastrow term (1.54) coincides with the solution of the two-body scat- 
tering problem (1.50) for a particular choice of variational parameters A,B,D. Notice 
that the parameters A, B, D are subject to a variational optimization within a quantum 
Monte Carlo framework, although the trial wave function, constructed from a two-body 
scattering solution, is fixed for any given choice of the mass of the particle and of the 
interaction strength. However in practice it might be convenient to keep the functional 
form of the solution and to treat its characteristic coefficients as variational parameters. 
A typical Jastrow term of the trial wave function of this symmetrized functional form is 
given in Fig. (2.1) in the Section, devoted to the Monte Carlo methodology. 

Worth noticing that in case of the Yukawa potential the pair distribution function at 
zero can be finite for r = 0, as it happens for the Coulomb potential. This can be easily 
confirmed if one takes a series expansion from Eq. (1.52). The typical value of the leading 
component A/D in the conditions of our problem is of the order of 10, therefore the zero 
value of the trial wave exp(— A/D) is very small and practically indistinguishable of zero, 
as one can see from Fig. (2.1). 

1.3.2.4 Scattering on repulsive van der Waals 

A similar derivation path can be applied to obtain the solution of the two-body scattering 
for the system of Rydberg atoms, interacting via the simple 1/r 6 interaction potential. The 
Schrodinger equation for the two-body Hamiltonian in the reduced units for this problem 
(see Sec. 5) in the system of the center of masses reads as 

— w"(r) H — qW(t) = Ew(r) (1-54) 

where we used a familiar substitution f(r) = w(r)/r. The finite low-energy solution of the 
latter equation has the following form (normalization is omitted): 



f(r) = ^=K -,^\ (1.55) 
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where K(a,x) stands for the modified Bessel's function of the second kind. Its expansion 
in the series of r is given by 



The short-range behavior of the function is dominated by the exponential term e^~2^~, 
which smoothly approaches zero as r — >■ 0. As in the case of the Yukawa potential, the 
interaction constant Cq and the power (-2) can be conveniently treated as variational 
parameters, while keeping the overall functional form of the trial wave function intact. A 
typical form of a trial wave function and a pair distribution function for the system is 
presented in Fig. (1.1). 
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Chapter 2 

Quantum Monte Carlo methods 



2.1 Introduction 

Quantum Monte Carlo methods (QMC) are very efficient and multi-purpose tools for the 
investigation of quantum many-body systems (for a detailed review of the methods see, for 
instance, [Cep95], [Gua98]). The use of Quantum Monte Carlo techniques provides a deep 
insight into the microscopic behavior of quantum states of matter. QMC are essentially 
ab initio methods, relying on a microscopic description of a system and gathering valuable 
information on its properties of the system through a numerical simulation. In certain 
cases, it turns out that this technique is the only tool for studying complex problems with 
reasonable calculation costs. In fact, in order to have a model, accessible for being solved 
analytically in the exact way, a physicist usually faces the necessity to make some kinds 
of approximations, but this can be avoided to a big extent by virtue of Quantum Monte 
Carlo simulations. As an example, the applicability of perturbation theory is limited by 
a small value of the perturbation parameter, while the QMC methods do not present 
any restrictions of this kind. Quantum Monte Carlo technique allow to find the ground 
state solution of the many-body Schrodinger equation at zero temperature. As follows 
from the name, the Quantum Monte Carlo methods are based on stochastic numerical 
algorithms of different sorts, that by nature are especially advantageous when the system 
in question possesses multiple degrees of freedom. As for any other stochastic procedure, 
the QMC methods provide results with a certain statistical error that can be diminished 
by performing longer measurement series. 

We are interested in studying the quantum properties of a given system. The quantum 
effects manifest the most when they are not disturbed by the thermal motion, that is at 
the lowest temperatures, when the system maintains in its ground state. For this case 
a possible method of choice to address the problem is the diffusion Monte Carlo (DMC) 
method. For a bosonic system this method allows to obtain the exact result for the energy 
of the ground state, as well as for any diagonal property of the state. 

In this chapter first of all we discuss the variational Monte Carlo method as the simplest 
one. Then we will discuss the bosonic diffusion Monte Carlo method and give a detailed 
explanation on how trial wave functions are constructed. Finally, we will present the main 
ideas for the implementation of the sampling of the quantities of interest. 
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2.2 Variational Monte Carlo 



2.2.1 General notes 

The simplest of the Quantum Monte Carlo methods is the variational method (VMC). 
The general idea behind this technique is to find an approximate wave function ipT, called 
trial wave function (or, sometimes, variational wave function), and then by sampling from 
the probability distribution 

p(r) = |Vw(r)| 2 (2.1) 

calculate averages of physical quantities. It can be shown that Et, the expectation value 
of the Hamiltonian, is an upper bound to the ground-state energy E . By expanding 
the normalized trial wave function in the basis of the normalized eigenfunction of the 
Hamiltonian 



oo 



V>t = I>^ ( 2 - 2 ) 

8=0 

oo 

EM 2 = 1 > (2-3) 

i=0 

one can rewrite the variational energy as 

Et - ii I / \ (2.4) 

oo oo oo 

= <£^|£|5>^> = E (2.5) 

i=0 j=0 i,j=0 

oo oo oo 

= E c*cAj E i = E E ih\ 2 > E ° E N 2 = E ° (2-6) 

i,j=0 i=0 i=0 

where i% stands for a corresponding eigenvalue (energy of the z-th state). By minimizing 
the variational energy with respect to the parameters entering into it one can optimize 
the wave function within the given class of considered wave functions. The only situation 
when the zero variance can be reached is when the wave function is exactly known. 



2.2.2 Usage of the VMC method 

If the wave function, corresponding to the ground state, is exactly known the sampling 
by means of the Variational Monte Carlo method permits to evaluate exactly any static 
property of the system within some statistical errors. Such systems are scarce; the ground 
state of a system of hard rods [KMJ99] and the Tonks- Girardeau gas [Gir60] are among 
the most famous ones. The Variational Monte Carlo method in this case can provide the 
correlation functions, which could be not accessible directly through the wave function. 
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The Variational Monte Carlo approach provides not only a valuable description of the 
quantum systems, but it also can be used as a first step to deliver a good quality input for 
the diffusion Monte Carlo method. The efficiency and even applicability of this method 
depends substantially on the optimization of the trial wave function within a chosen class 
of functions. 



2.2.3 Notes on algorithmic realization 

Let us stick to a coordinate representation in the following description, since it is the easiest 
way to work with external or interparticle potentials. Consider a common D— dimensional 
Euclidian space with a system of N p particles inside. The probability distribution function 
in this system will be a function of N ■ D variables p(R) = rjv_). The mean value 

for an arbitrary operator Z is therefore calculated as an integral of N ■ D dimensions in 
the following form: 

, z s = [ ■•■ J Z(t u rjy p )p(ri, r Np )dr 1 ...dr Np 
I ...fp(r 1 ,...,r Np )dr 1 ...dr Np 

It is clear that the complexity of the estimation of this integral with conventional non- 
stochastic methods, based on a grid calculation, grows extremely fast with the number of 
particles, and already for a few dozens of particles its calculation is unreachable. On the 
other hand, the stochastic procedure on which the Monte Carlo methods are based is not 
affected strongly by the growth of the dimensionality of the problem. The basic idea of the 
variational technique is to generate configuration states Ri,...,R m with the probability 
distribution p(R) and by means of these states estimate the average value of the operator 
(Z), 

i m 

<Z>«-£2(H*). (2.8) 

Every state is obtained only from its preceding configuration, thus this set of states is 
indeed a Markovian chain. The Metropolis algorithm [MRR + 53] despite its simplicity is a 
very efficient tool to produce the chain with the desired probability p(R). The transition 
from the old R configuration to the new one R' is accepted with the probability P(R — >• 
R'), expressed by the formula 

P(R ^ R) = \ 1, ifp(R')>p(R) (2 ' 9) 

In a quantum system, the probability distribution is given by the square of the wave 
function. The way we construct the particular trial wave function for different quantum 
systems will be discussed in Sec. 2.4. 

There are distinct approaches to perform transitions (or trial moves) between different 
particle configurations. A straightforward way for creating a new state is to move one 
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particle at once or all the particles at once r^ = r.; + pi, i = 1, N p , where the displacements 
Pi are a set randomly chosen vectors with a certain upper bound < U. The limiting 
amplitude U can be adjusted in order to have a desired acceptance rate. It is readily seen 
that for very small values of U the trial moves are accepted almost all the time, but the 
adjacent states are strongly correlated, that affects the variance of the overall statistics, 
and in the limit U — > all the states are the same, making the whole calculation pointless. 
On the other hand, the steps of very large amplitude U are accepted only a small fraction of 
time, which also leads to a poor performance. The value of U can be optimized to ensure 
the maximum total displacement of the whole system by means of a set of benchmark 
calculations, but in practice a simple rule of thumb to have the acceptance rate of the trial 
moves in the range 0.5 — 0.65 provides good enough results. 

A generalization of these two strategies of particle displacements, when only a particular 
fraction of the points is shifted together, can bring even faster performance. The group of 
particles to move in one step can be chosen randomly. The advantage of this technique lies 
in the possibility of fine tuning the calculation parameters in order to achieve an optimal 
point in terms of the interplay between the correlation of the states and the acceptance 
rate. 



2.3 Diffusion Monte Carlo technique 

The diffusion Monte Carlo method (DMC) is a stochastic computational technique applied 
to systems at zero temperature, when all of the thermal motion can be neglected. The key 
point of the DMC method is to provide the solution of the time-dependent Schrodinger 
equation in imaginary time, which is known to exponentially decay to the ground state 
solution in the limit of long times. By means of the diffusion Monte Carlo method the 
equation of state for the system as well as diagonal properties can be calculated exactly 
with the only cost of controlled statistical noise. 



2.3.1 Schrodinger equation 

The wave function of a quantum system obeys the Schrodinger equation 

ih^(R,t) = HxP(R,t) . (2.10) 

Our aim is to find the ground state properties of the system, rather than its actual time 
evolution. By substituting the time variable by an imaginary one r = —it/h, one arrives 
to another representation of the Schrodinger equation 

- J%(R, t) — (H — E)i/;(R, r), (2.11) 

where E stays for a constant energy shift close to expected ground state energy. The 
latter equation (2.11) has a formal solution ip(R,r) = ' t i/j(R,0). One can expand 
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this solution as the sum over the eigenstates of the Hamiltonian H(j) n = E n (j) n , with the 
eigenstates taken such that the eigenvalues are growing with ascending indexes, that is Eq 
is the lowest among the eigenvalues. Performing the expansion, 

oo 

V>(R,r) =5><MR)e- (i ^ )T . (2.12) 

i=0 

It can be easily seen that the exponents in the sums behave differently (decay or grow) 
provided the sign of (Ei — E) is positive or negative. For large enough r the only component 
of the sum (2.12) that survives is the one, corresponding to the ground state. All the other 
terms decay in time exponentially fast (we suppose the spectrum to be discreet): 

^(R,r) -»• a 0o (R) e- T{E °- E \ if r ^ oo (2.13) 

A general expression for the Hamiltonian of a many-body system of N p particles, sub- 
jected to an external force field V ext {ri), depending on a particle's position, and interacting 
internally through a pairwise potential VJnt(|i"i — r^j) can be given as 

fc2 N p N p N p 

H = ~ E A? + E + E VWk< - r il) • (2-14) 

/m %=i i=i i<j 

We take a notation for the constant D = h 2 /2m, referring to D as a "diffusion" constant, 
the sense of which will be clarified later. The Schrodinger equation (2.11) reads as 

- ^(R, r) = -DA R ^(R, r) + ( WR) " E)4>(R, r), (2.15) 
or 

with the label R in the Laplacian denoting the differentiation with respect to each scalar 
component of the vector R. For shortness, the summation over the internal and the 
external potentials is denominated by a single term Vtotai(R) = Kxt(r^) + ^int(| r iil)- 

The last term in Eq. (2.15) (V totau \(R) — E)ip(R, r) is diagonal, and it affects the normal- 
ization through a specific step, called branching. As we mentioned above, a key ingredient 
of the diffusion Monte Carlo technique is the use of a trial wave function which allows to 
radically reduce the calculation efforts needed to reach a required accuracy in the result. 
Hence, one gives an approximation of the true ground state solution ipo(R, r) by a certain 
trial wave function ^t(R, t), which is subject to a correction inside the algorithm by means 
of the branching. The whole approach, referred to as importance sampling, stems from an 
analog of the Schrodinger equation for the product of the true wave function and the trial 
wave function 

/(R,r) = ^(R)^(R,r). (2.16) 

The expectation value of an operator Z over this product of the wave functions can be 
thought of as a mixed estimator (Z) = J ^Zip dR/ / ipr^ dR. Performing a substitution 
of / inside the Schrodinger equation (2.15) the following equality can be drawn 

- J^/(R, r) = -DA R /(R, r) + (£ local (R) - E)f(R, r) + DV R (/(R, r)F(R, r)) (2.17) 
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where E local is the so-called local energy, that is the expectation value of the Hamiltonian 
with respect to the trial wave function (a broad discussion on the calculation of the values 
is given in Section 2.5). 

E^(R) = = ( W R) - E) - D^S*^ . (2.18) 

The action of the Metropolis algorithm in the variational Monte Carlo method, that is 
the averaging of an operator over the probability distribution p(R) = |^t(R)| 2 ; is therefore 
the averaging of the local energy with respect to the set of N c particle configurations 

i N c 

E = — ££ local (R,). (2.19) 

Notice that the local energy is equal to the eigenvalue of the Hamiltonian, when the trial 
wave function coincides with the eigenf unction. It means that E local is a constant ground 
state energy with zero variance if "0t(R, t) = ^>o(Rj r )- The deviations of the trial wave 
function from the exact solution result in a growth of variance in the local energy, which 
can serve as a quality criterion for the trial wave function. The notation F stands for 
the drift force, a vector value, equal to the gradient of the field defined by the trial wave 
function, multiplied by a convenient factor, 

2V r V>t(R) 

F = — RV Z\ ; . 2.20 

The drift force moves the walkers {Ri} towards the region where ipT is large. 

Notice that the probability distribution function in a classical system is defined by a 
Boltzmann factor /(R) ~ exp(— £ , pot (R)/A; jB T), with E pot (R) standing for the potential 
energy in the system. The force, produced by the field, is equal to its gradient with the 
negative sign, that is F(R) = — VE pot = Vlnp(R). Making a formal substitution of 
the trial wave function instead of the probability density p(R) = ?/4(R), we can recover 
precisely the same form of the force, as in the definition (2.20). 



2.3.2 The Green's function 

The Schrodinger equation (2.15) can be resolved formally in the following way: 

(R|/(r)> = £(R|e-^)|R>(R|/(0)>, (2.21) 

R' 

where / stands for time-dependent wave function appearing in the equation. 
The term (H\e~^ H ~ E ^ T \Ii.') is the Green's function G(R, R',r) of the Hamiltonian, acting 
as a propagator of the system. The Schrodinger equation can be rewritten in the integral 
form in terms of G(R, R', r) as 

/(R, r) = J G(R, R, r)/(R, 0) dR (2.22) 
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The last integral equation is an analog of the traditional form of the Schrodinger equa- 
tion (2.10) but it allows to obtain the solution /(R, r) by virtue of a multidimensional 
numerical integration, where the Monte Carlo methods are applicable. The exact form of 
the Green's function is not known, however it can be expanded in power series of r, there- 
fore with a certain approximation the solution can be reached in a number of integrations, 
depending on a the time-step size. Hence, the accuracy and the overall computational 
complexity to obtain the final solution depend on the time step. The solution of Eq. (2.22) 
after a single time step is then given by 

/(R, r + At) = J G(R, R, Ar)/(R, r) dR (2.23) 

As it was stated before, after a large enough number M of time steps the solution 
of Eq. (2.22) decays to the ground state, while all the other components of the solution 
exponentially disappear. 

/(R, MAr) -»• ^t(R)0o(R), when M oo (2.24) 

It is natural to separate three different operators inside the Hamiltonian 

H = (—DA) + D((V R F) + FV R )) + (-E + £ local (R)), (2.25) 

that is 

H = H x + H 2 + H 3 , (2.26) 

and write down the Green's function for each of the operators Hi, 

Gi(R,R',r) = (R|e- T ^|R) . (2.27) 

If components of an operator do not commute, its Green's function (or exponent of the 
operator) cannot be represented as a product of the Green's functions of the components. 
This is evidently the case for the components Hi of the Hamiltonian. Nonetheless, there is 
a possibility to use approximations of different orders of r for exp(— tH). The first-order 
approximation for the exponential of H corresponds to 

e- rA = e -^i e -^ e -T# 3 + ( r 2) _ ( 2>28 ) 

The integration of this formula, when the term 0(t 2 ) is neglected, yields the following 
expression for the Green's function 

G(R, R, r)= JJ Gi(R, R 1; r)G 2 (R 1 , R 2 , r)G 3 (R 2 , R', r) dRi dR 2 (2.29) 
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Notice that in accordance with (2.27) the Green's function G is a solution of the Bloch 
equation 

-j^G(R,R' J T) = ^G(R,R,r), i = 1,2,3 (23Q) 
G(R,R,0) = (J(R-R') 

The first (kinetic) component G\ satisfies the equation 

<9Gi(R,R',r) 



dr 



£)AGi(R,R',r) (2.31) 



which coincides formally with the diffusion equation. The diffusion constant of this equa- 
tion is equal to D, and this is the reasoning of the name that we adopted for the constant. 
The last equation (2.31) can be easily solved in the momentum representation, as the ki- 
netic energy operator is diagonal in this representation. Rewriting the solution again in 
the space coordinates, we get a well-known formula 

G 1 (R, R, r) = - 1 e-* 3 ^ (2.32) 

The second term G2, containing the drift force, satisfies the equation 

0G 2 (RJt',r) = jDVr(G2(R) r , r)F) {2 33) 

is also easily resolved 

G 2 (R,R,r) = 5(R-R c (r)), (2.34) 
where R c (t) stands for the solution of the equation 

; ^ = df(rjt)), (2 35) 

R o (0) = R' 

which defines the motion of the system, subjected to the drift force F. 
The third equation from the set (2.30) is evident to solve: 

G 3 (R, R, r) = exp{(£ - £ local (R)) r} 8(R - R) (2.36) 



which is generally referred to as the branching component, as it controls multiplication 
and annihilation of walkers in the algorithm. 
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2.3.3 The primitive (first-order) algorithm 

If the wave function of the system ip(R, r) is real and positive, as it happens for the ground 
state of a Bose system, it can be treated as a population density distribution: 

-^walkers 

V(R,r) = C £ 5(R-Ri(<r)), (2.37) 

i=l 

with a constant C > and Ri(r) are time-dependent positions of a single particle set, 
referred to as a walker. The formula (2.37) should be taken in a statistical sense, the average 
of any value A over the left hand side and right hand side distributions are equal to each 
other in the limit when size of the population Nw tends to infinity / A(R)/(R, r) dR = 

N w 

lim / i4(R) C6(R — Rj(r)) dR. The walker resides in the coordinate space of 3N 

iV u ,->oo j = i 

dimensions, and the infinitesimal probability /(R, r) dR is equal to the probability to 
encounter a walker in the infinitesimal range dR near the point R in the moment r. 

The algorithmic implementation of the time evolution of the system according to its 
Hamiltonian (2.25) is an evolution of a system of walkers upon the action of the elementary 
components of the Green's function, that is the transition matrices (2.32), (2.34) and (2.36). 

The first function Gi determines a diffusion of the whole system of walkers according 
to the equation 

Ri(r + At) = a + R(r), (2.38) 

where a stands for a displacement, taken randomly from the three-dimensional normal 
distribution exp{— a 2 /(4ArD)}. 

The second function Gi gives rise to the drift of the set of walkers towards areas in the 
configuration space where the trial wave function is large: 

R 2 (r + At) = R(r) + (DAt)F(R) . (2.39) 

Notice that the Green's functions of the first two kinds (2. 32), (2. 34) have a normaliza- 
tion / G 12 (R, R', t) dR = 6(Rf), that is on each of the steps the number of walkers stays 
unchanged. 

This is not the case for the third propagator, corresponding to the branching term, 

/( 3 )(R, r + Ar) = /(R, T ) e -(E-E^\R)) Ar ? {2AQ) 

since the integral over its Green's function G 3 (R, R',r) (2.36) is clearly not equal to 
unity. The most evident interpretation of the action of the third propagator is that each 
walker has a certain value attached to it, commonly thought as its "weight". This value is 
recalculated on each step for every walker. The formula (2.40) suggests that the walkers 
of less _E local are favored and the contrary are disfavored. A clear disadvantage of this 
scheme is that one would wish to have a better statistical representation for the favorable 
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areas of the configuration space, however the generation and elimination of walkers is not 
possible. The alternative and widely applied way to overcome the problem is that the value 
c = exp{-(E - £ local (R)) At} is treated number of exact copies of the walker. As 
one can see, this is a non-integer number that cannot be readily used for it. However, one 
can apply a randomization by throwing a random number from the uniform distribution 
[0, 1) and afterwards take [c] + 1 duplicates of the walker, if the number is within [0, {c}], 
and [c] is otherwise (the standard bracket notations [.] and {.} stand for the integer and 
the fractional parts, respectively). The relation ([c] + l){c} + [c](l — {c}) = c ensures the 
correct mean action of the propagator. 

The equation (2.40) also suggests that the branching can be effectively controlled by a 
dynamic choice of the reference energy E, that can be essential to avoid a collapse of the 
set of walkers as well as its undesirable expansion, that can make the simulation stall. 

Looking closely one can notice that the diffusion and the drift steps represent only the 
sampling from the trial wave function. It means that being applied without the branching 
term, they are equivalent to the application of the variational Monte Carlo method (see 
[HJR94] and Section 2.2). The third branching step ensures that the system of walkers 
"prefers" the areas of higher ipo (it is often called a correction of the trial wave function), 
and the overall sampling is taken from i^t^o rather than As it was commented above, 
if the trial wave function coincides with the ground-state solution (or, more generally, with 
an eigenfunction of the Hamiltonian) , the local energy becomes equal to the eigenvalue, 
hence the branching factors are the same for all the walkers. The action of the branching 
step in this case does not affect the final result. 



2.3.4 Second-order algorithm 

In the previous section, we have described the simplest (first-order in r) approximation 
(2.28) of the Green's function. The order of the application of the three propagators is 
clearly irrelevant, since overall results (for instance, the energy of the ground state) will 
depend linearly on the time step r. To obtain the final result for a quantity of interest 
one should take a short enough r to move the time-related error below the statistical 
noise or perform a series of simulations with distinct time steps, and then find a linearly 
extrapolated value. In practice the second approach is much more practical, since the 
time bias of the result can be very pronounced. However, a clear drawback is that the 
correctness of the linear dependence might be valid for undesirably small times. This is 
where the higher-order algorithms become useful. 

The second-order in r expansion of the exponent in the Green's function can be given 
in the form 

which is not the only possible way of representation, but probably the most efficient for 
actual application [Chi90]. The final result for the equation of state of the system, yielded 
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by the formula (2.41), does not have linear in r elements, thus for small enough time steps 
the dependence is quadratic. In this case one can again perform an extrapolation to r — > 
via a series of measurements, or alternatively, find a time step short enough to guarantee 
the smallness of the time step-related error in accordance with a required accuracy level of 
the simulation. This accuracy is generally described by the statistical variance of the result. 
The choice between the two approaches must be taken regarding the interplay between the 
additional calculation costs, related to the extrapolation, and the slower evolution of the 
system in case of a shorter time step. 

Let us explain the second-order algorithm, that relies on the expansion (2.41). Each 
propagation in time corresponds to a constant shift At/2 with respect to the current time. 
The state of a displaced walker is changed on every step from TU-i to Rj. The effect of 
the branching propagator is a corresponding multiplication or suppression of a walker in 
question. Since the walker is moved in a loop, the choice of the first step is arbitrary, and 
the list of the operations can by reordered. 
The calculation procedure: 

1) The first propagator, random Gaussian move (2.32): 
R = Rj_i + a , a is taken from exp(— x 2 /At) 

2) The second propagator, drift with F(R) (2.34): 
R = R + F(R) Ar/2 

R" = R + (F(R) + F(R) At/4 
R'" = R + F(R") At 

3) The branching propagator (2.36): 
R'" unchanged. 

2.4 Constructing the trial wave functions 

2.4.1 Motivation 

In this Section, the development of trial wave functions for bosonic systems is discussed 
from the technical point of view. The purpose of the following discussion is to present the 
theoretical basis of the construction of Nosanow-Jastrow trial wave function of a general 
kind in liquid and solid phases, its relation to the actual implementation, the technical 
issues that one faces, and their solutions. 

2.4.2 Nosanow-Jastrow trial wave function 

The natural way of constructing a bosonic trial wave function is to take it in the form of 
a product of one- and two-body correlation terms: 




(2.42) 



j<k 
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This form of the trial wave function is generally referred to as the Jastrow (or Bijl- 
Jastrow) trial wave function, first proposed by Bijl [Bij40] and then by Jastrow [Jas55]. The 
pairwise interaction of particles in the system is taken into account by the two-body Jastrow 
term /2( r )- The pair correlation is clearly lost for large enough interparticle distances, that 
is this term must approach 1 asymptotically. The use of the periodic boundary conditions 
puts additional restrictions to fzi?). Namely, in order to avoid additional contributions 
to the kinetic energy from distances r > d p \ anc (rfpiane stands for a distance to the plane, 
closest to the center of the simulation cell), the Jastrow term should be constant in that 
region. For a rectangular simulation box it is one half of a minimum box dimension 
^pianc = min(L :r , L y , L z )/2. For a calculation cell of the shape of truncated octahedron (see 
Appendix B) this distance is equal to 

rf planc = ; : === . (2.43) 

"I - 1 IT2 "T" 



l/Ll t i /L 2 -r 1/L 2 

The value of the Jastrow-Bijl term at distances larger than the cut-off one r > L/2 is 
not necessarily f 2 (L/2) = 1, but, in principle, can be any arbitrary constant /2(L/2) = 
const, r > L/2. Still, the unitary value might be convenient to use in the evaluation of 
the product of Bijl- Jastrow terms n /2( r i — r j), which in code is usually implemented as 

i<j 

evaluation of exp m /2( r « ~ r j)- With the choice / 2 (L/2) = 1, r > L/2 the contributions 

i<j 

of pairs with |r.j — r ^ | > L/2 to the sum is zero. 

The one-body Jastrow term /i(r) is introduced to take into account an external po- 
tential, present in the system. It can also define symmetry properties of the system, for 
example, the localization of the particles in crystalline nodes. The form of this one-body 
term is typically obtained from a solution of the Schrodinger equation for a single particle 
in the chosen potential. 

For a quantum Monte Carlo simulation of a solid phase one might need to induce a 
corresponding crystalline symmetry to the trial wave function. This is done by "pegging" 
the particles to the nodes of a crystal through a multiplication of the Jastrow term by a 
particular factor, depending on the coordinates of the particles. A straightforward way to 
realize such a condition is to consider a one-body term 

/i(r i )=^(|r i -6|) (2-44) 

with the configuration {£} standing for a set of crystal nodes' positions and g(r) is a 
function, which localizes each particle to the site £j. This factor is generally referred to 
as Nosanow term 1 , and a corresponding trial wave with a Jastrow two-body term 

* NJ (r 1; r Np ) = n f2(r ij ) J] 9(1* - 61) (2.45) 

i<j k=l 



1 A not so brief explanation of the efficient construction of crystalline guiding wave functions can be 
found, for instance, in the recent work of Cazorla and collaborators [CACB09]. 
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is called a Nosanow-Jastrow (NJ) trial wave function. The localization factor g(r) is usually 
chosen as a Gaussian with a localization parameter 7 = l/(2cr 2 ): 



since the Gaussian is a quantum mechanical solution for a 3D harmonic oscillator and 
usually is a good guess for describing a wave function close to a potential minimum. The 
parameter 7 is determined through a variational optimization in VMC calculations. 

Quantum Mechanics requires that the bosonic wave function is symmetric with respect 
to an exchange of equivalent particles, and that the Nosanow term does not satisfy this. 
Evidently, the analytical form of the symmetric trial wave function in this case becomes 
quite cumbersome for implementation and extremely expensive in terms of permanent cal- 
culation time, as it should in general contain a sum over all the permutations of the lattice 
nodes in the system. Nonetheless, the results of the quantum Monte Carlo simulations 
when the Nosanow-Jastrow term is used are for practical purposes indistinguishable from 
that of the symmetrized trial wave function [CB08b, CACB09], as far as the energetic 
properties of a crystal are concerned, and the exchange energy is usually negligible. In 
this Thesis we will use this form of the NJ trial wave function throughout all the quantum 
Monte Carlo calculations of solid phases. However, the physical quantities, related to Bose 
statistics, may not be treated with the NJ trial wave function, since any particle exchange 
is suppressed. 

The quantum Monte Carlo technique requires knowledge of not only the trial wave 
function itself, but also of its first two derivatives. In the case of the DMC method the 
second derivative enters in the calculation of energy, which is needed for the branching term. 
Hence the second derivative of the trial wave function is required not also for evaluating 
averages of the observables, but even for the time evolution of the system. 

The actual implementation can be substantially improved, if one takes into account 
that the trial wave function appears in the implementation of the method in three distinct 
combinations, which can be presented as functions and can be tabulated. 

A. The logarithm of the Jastrow term (needed in Metropolis algorithm in VMC simula- 
tions and also for estimations of the non-local quantities, for instance the one-body density 
matrix) 



B. The logarithmic derivative of the trial wave function, which is required in the cal- 
culation of the drift force (2.20)) 



0(1**- 61) = e 



(2.46) 




(2.47) 





C. The second derivative of the Jastrow term appears only in a linear combination with 
the drift force in the calculation of the kinetic energy. The following representation takes 
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place: 



,io„i,_> /"M , /7'MV , ™K,.M d-i/'W 



*~m HtMJ + t - — tM- (2 ' 49) 



with d stands for for a dimensionality of the problem. 



2.4.3 Explicit expressions for wave functions 
2.4.3.1 Trial wave function for hard sphere potential 

The problem of scattering on a simple hard sphere potential, Yukawa potential and common 
potential between Rydberg atoms was discussed in Sec. 1.3.2.2. The hard sphere and the 
considered model potential of Rydberg atoms make the wave function vanish when two 
particles meet each other, which means that three-body collisions are greatly suppressed. 
The same is true for a Yukawa potential at low density due to its similarity to the Coulomb 
potential. If three-body correlations are neglected, at small interparticle distances the 
two-body Jastrow term /2(r) can be conveniently approximated by the solution f(r) of 
the two-body scattering problem, that is by the wave function of a system of two particles. 
At large distances, the pair wave correlation function asymptotically tends to a constant, 
as the particles lose correlation. 

Taking these facts into account we introduce the trial function in the following way [GBC99] 
(here we use a dimensionless notation in which the distances r are mesured in units of the 
hard sphere radius cl^d and the energy E is measured in units of h 2 j \ma\ D )) 

Asin(V2E(r - 1)) 



I? exp \ \ , \r\> R 



in 



The Jastrow term has to be smooth at the matching point R m , that is 
A. the function /2( r ) itself must be continuous: 



Rm ^ ^ 

B. derivative f^ir) must be continuous 

A^2E cos{V2E{R m - 1)) Asixi(yfE(R m - 1)) B ( R 



^exp(-^l (2.52) 
a la) 



C. the local energy f 2 (r)~ 1 (—h 2 Ai/2m — h 2 Ai/2m + Vj nt (rj — rj))/ 2 (r) must be contin- 
uous 

Rr 



B exp 

R_ 

a 



' in. 



2E = Rm(x) - p \ a 7 (2.53) 

1 — B exp 
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The solution of this system is 

r e-^ 

~ sm(u(l-l/R))e-2^ + u^ (2 54) 

B _ u 2 exp(0 



where we used the notation u = \/2ER and £ = R/a. The value of £ is obtained from the 
equation 

11 u(£-2) , 

1 - — = - arctan V ^ - ' 2.55 

There are three conditions for the determination of five unknown parameters, conse- 
quently two parameters are left free. The usual way to define them is to minimize the 
variational energy in variational Monte Carlo which yields an optimized trial wave func- 
tion. 

2.4.3.2 Trial wave function for Yukawa potential 

The construction of a trial wave function for a Yukawa system can be done in different 
ways. The first one, widely employed in our calculations, is a use of the solution of a 
two-body scattering problem, that yields a satisfactory short-range approximation for the 
trial wave function, which is valid for a dilute system. At large distances the trial wave 
function is intended to arrive smoothly to unity at the half size of the simulation box 
(see Sec. (1.3)). If the Jastrow term is chosen in a form /(r) = e~ u ^ [RC67], this can 
be achieved by a symmetrization of the trial wave function with respect to the inversion 
r — > (L — r) as 

u(r) := u x {r) + m(L - r) - u x {L/2) (2.56) 

that brings the logarithm of the Jastrow term and its first derivative in the point L/2 to 
zero. 

According to Eq. (1.51), the solution of the two-body scattering problems reads as 

, x , ( h (const \fr)\ i—r A A 2 „ „. 

«i (r) = In f u ^ 1 = In v 7 ! + -r - — r 2 + C(r 3 ) (2.57) 

where A is a constant, subject to optimization. This solution formally coincides with the 
scattering solution provided A = Vq (Vq stands for the interaction strength constant of the 
original two-body scattering problem (1.49)). An optimal value of A should therefore be 
close to Vo. 

The most productive way to generate trial wave functions in the case of the Yukawa 
potential appeared to be the hypernetted chain method [CP79], based on an iterative 
solution of a set of Euler-Lagrange equations 

« Wflf -M-l,-* (2,8) 
du n (ri, ...,r„) 
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E/N 


VMC Jastrow TWF 


19.862(3) 


VMC HNC TWF 


19.634(3) 


DMC result 


19.595(3) 



Table 2.1: The energy per particle without size correction for different t. w. functions 
compared to the exact DMC result. 



The more detailed explanation of the method is given in Section 2.4.3.4. 

A comparison of two-body correlation factors from Eq. (2.57) and the HNC solution is 
given in Fig. (2.1). The particular conditions of the data correspond to VMC and DMC 
calculations with 64 particles in the truncated octahedron cell, A = 0.46, p = 0.024 (for 
details on the used model and involved parameters see Section 4). 




r 

Figure 2.1: Two-body correlation factor for a system, interacting via a repulsive Yukawa 
potential in the liquid phase, in the Jastrow parameter-optimized model (red solid line) 
and the HNC one (green dashed line). The particular conditions of the data correspond to 
a DMC calculation with 64 particles in the truncated octahedron cell, A = 0.46, p = 0.024. 

HNC f(r) leads to better estimates of the energy; a simple comparison is shown in 
Table (2.1) (the conditions of the simulation are the same as in the figures above, no finite 
size correction added). 



2.4. Constructing the trial wave functions 



39 



2.4.3.3 Trial wave function for repulsive van der Waals 



In our simulations of the bulk system with the pairwise van der Waals interaction l/|r| at 
zero temperature we use the short-range approximation of the two-body scattering prob- 
lem, as it was solved in Sec. (1.3). The technical procedure to obtain the functional forms 
of the solution follows the derivation, used in the case of the Yukawa potential. The final 
result for the Jastrow factor reads as a logarithm of the first term in the expansion (1.55): 

Mr) = V" ln ^-— (2 ' 59) 

u(r) = ui(r)+u l (L-r)-2u l (L/2) . (2.60) 

In this equation the factor C§ comes from the interaction strength and is constant. 
Nonetheless, it can be treated as a parameter and variated in order to optimize the trial 
wave function by minimizing the VMC energy. Typical forms of the guiding wave functions 
and pair distributions are presented in Fig. 2.2. 
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Figure 2.2: A typical two-body correlation factor and pair distribution function for a 
system, interacting via a repulsive van der Waals 1/r 6 potential in the liquid phase. The 
particular conditions of the data correspond to a DMC calculation with 108 particles in 
the cubic cell, p = 6.7. 



A similar result is provided by the cusp condition demand, when the leading term of 
the guiding wave function is sought in a suitable exponential form f(r) = exp(— a/r 6 ), 



40 



Chapter 2. Quantum Monte Carlo methods 



with parameters a, b. After a substitution of this functional form into Eq. (1.54)) one finds 
that the equation can be satisfied only when b = 2, while a is still arbitrary. It is easily 
seen, that this procedure yields the leading term of the two-body scattering solution. 

2.4.3.4 Hypernetted chain method 

The hypernetted chain method is a technique to solve many-body problems in homogeneous 
and inhomogeneous media [Kro98],[vGd59]. In this scheme, the static structure factor S(k) 
that minimizes the variational energy in the subspace of Jastrow wave functions has the 
form 

S(k) = m , (2.61) 

\Jt 2 {k) + 2t(k)V ph (k) 

with t(k) = H 2 k 2 /2m and V p h(k) the so-called particle-hole interaction. Its Fourier trans- 
form FT[V p h(k)] = V p h(r) satisfies the following equation in coordinate space 

V ph (r) =g{r)V{r) + ^-\V^{r)\ 2 +{g{r)-l)uj I {r) , (2.62) 

where V(r) and g(r) are the bare two-body potential and the pair distribution function 
(the Fourier transform of S(k)), respectively. Finally, in momentum space the induced 
interaction uj(k) becomes 

WlW = _^ W IWM. (2 , 3 ) 

In this way, Eqs. (2.61), (2.62) and (2.63) form a set of nonlinear coupled equations that 
have to be solved iteratively. The Fourier transform of the resulting S(k) provides g(r) 
and, in this scheme, the optimal two-body Jastrow factor results from the corresponding 
HNC/0 equation 

/ 2 (r) = JgV)e- N(r)/2 , (2.64) 



where N(r) is the sum of nodal diagrams, related to S(k) in momentum space by the 
expression N(k) = (S(k) - lf/S(k). 

2.5 Estimators for physical quantities 
2.5.1 Local energy 

The general form of a Hamiltonian of a system of N interacting bosons in an external field 
is (2.14): 

t2 N N N 

H = + E^xt(r) + E^ntdr, - r fc |), (2.65) 

i=l i=l j<k 
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where m is mass of a particle, V ex t(r) is the external field, V| n t(|r|) is the two-body particle- 
particle interaction potential. Given the many-body wave function \I/(ri, ...,rjv) the local 
energy is defined according to (2.18): 



E loc \ 



ri, ...,r N ) 



H^(r 1 ,...,r N ) 



(2.66) 



The external field and particle-particle interaction are diagonal in this representation 
and are calculated trivially as a summation over particles and pairs of the second and third 
terms of (2.65). Calculation of the kinetic energy, first term of (2.65) is more involved, as 
the Laplacian operator is not diagonal. 



2.5.1.1 Local kinetic energy 

In this Section, we will find the expression of the local kinetic energy 

h 2 A\&(ri, r N ) 



T loc (r u ...,r N ] 



2m tf(n,...,rjv) 



(2.67) 



Let us calculate the second derivative in two steps, as the first derivative is important 
for the calculation of the drift force. We consider the Jastrow form (2.42) of the trial wave 
function and will express the final results in terms of one- and two-body terms fi and fi- 
The gradient of the many-body trial wave function is given by 



Vi#(ri,...,r 



Nj 



#(ri,...,r 



N 



/lfo) r 

The full expression for the Laplacian is 

' flirt) v 



f[(r t ) r i f' 2 ^ Vi ~ rfc D Ti ~ rfe 



i /2(ki - r fc |) |r< - r fc | 



(2.68) 



A r .*(n,...,r 



*(ri,..., rjv ) 
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*(ri,...,rjv) 



| y ^(l r » ~ r fcl) r i - r k 



+ 



N 

E 



/2(|r* — rfcl) \/2(|ri-rjfe| 



(2.69) 



The kinetic energy can be written in a compact form 



T loc (r l7 ...,r w ) 



fc2 I N jv 

E + 2 E 4° C (|r, -r fe |)--E r N ) 

[i=l j<k i=l 



N 



N 



(2.70) 



42 



Chapter 2. Quantum Monte Carlo methods 



where we introduced notation for the one- and two-body contribution to the local energy 
(see, also, (2.49)) 

£! ° c(r) = -§M§!) 2 (27i) 

and introduced the drift force (see (2.48)) 

^(rx, r N ) = 2 ( + £ ffj* ~ r f r * ~ **) (2.73) 
\fi(ri)ri %£M\Ti -r k \) |r< -r fc |y 

2.5.1.2 Exponentiation 

It is convenient (see Eq.2.47) to do the exponentiation of the one- and two- body terms 
«i(r) = ln/i(r), u^ir) = ln/ 2 (r). The point is that numerically a better precision is 
achieved by working with numbers of the same order. The formula for the kinetic energy 
becomes more simple 

h 2 f N N 1 N 1 

T loc (r 1; ...,r N ) = -—\J2 u?(r,) + 2 £ <(|r, - r fc |) + - £ |F*(r lf r^)| 2 (2.74) 
zm [i=i j<k 4 i=i J 

with 

TV 

Fi(n, r,v) = 2ui( ri )^ + E«2(l r i " ffcl) ,^"^, (2.75) 
2.5.2 Static structure factor 

The static structure factor S(k) is the correlation function of the momentum distribution 
between elements —k and k (1.28): 

NS(k) = (p_ kPk )-\(p k )\ 2 , (2.76) 

Using the properties of the Fourier component p_ k = (p k )* it can be rewritten in a 
different way 

NS(k) = (\ Pk \ 2 )-\(p k )\ 2 . (2.77) 

The density distribution in coordinate space is the sum of 5-functions located at the 
positions of the particles: 

1 N 

n(r) = v Y^S(r-r l ). (2.78) 

i=l 
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By means of its momentum space representation (i.e. its Fourier transform) 



N 



N 



N 



Pk = Y. e ^ = H cos kri-iY, sin kr t 



i=l 



i=l 



i=l 



we obtain a simple expression for the static structure factor 



NS(k) = ({fl c °s kr^j + (j2 sin kr^j 



N 



cos kr i 



\i=l 



N 



sin kr j 



u=l 



In a homogeneous system the two last terms in (2.81) are vanishing, that is 



NS(k) = / fe cos kr^j + (j2 sin fcr * 



(2.79) 



. (2.80) 



(2.81) 



If periodic boundary conditions are used, the values of momenta are quantized and 
depend on the size of the box 

(2.82) 



x,y,z ■ 



2.5.3 Calculation of one-body density matrix 

The one-body density matrix (OBDM) gi of a homogeneous system described by the 
many-body wave function ip(ri, r^) according to (1.18) is equal to 



2i(|r'-r"| 



N 



J ... / tp*(r ', r 2 , r N )tp(r ", r 2 , r N ) dr 2 ...dr N 



J ... J |*0(n, ...,rjv)| 2 drx.-dr 



(2.83) 



N 



Since in DMC one does not sample directly the ground-state probability distribution 
0o, but instead the mixed probability ipT<Po (2.16) one obtains the mixed one-body density 
matrix as the output 



. J ... J ipUri + r, r 2 , rjv)0 o (ri, r 2 , r^v) dr 2 ...dr N 



J ...J ^{r 1} •••) r iv)^o(ri, -, rjv) dr x ...dr N 
This formula can be rewritten in a way convenient for the Monte Carlo sampling: 



(2.84) 



mixed / \ / - /[#(gi + r > r 2, r jv )(^(r 1 , r 2 , rjy)) 1 ]/(r 1 ,...,r JV )dri...dr JV 

yi l r J _ n 



S ...J f(r 1 ,...,r N )dr 1 ...dr 



N 



(2.85) 



where we have used the asymptotic formula (2.24) and have taken into account that in a 
homogeneous system g q depends only on the module of the relative distance. If the trial 
wave function is chosen as a product of pair functions (2.42) then using the notation (2.47) 
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w(|r; - Tj\) = ln/ 2 (|ri - Tj\)) and f\ = 1 one has i/; T (ri, ...,r N ) = U exp{u(\ri - r^)}. 
Then the ratio of the trial wave function appearing in (2.85) becomes 

= e 3 tp(s;ti(|r 1 + r-r i |)-ti(|r 1 -p i |)l (2.86) 



■4> T (ri, -,r N ) 



J> 1 



In order to gain better statistics it is convenient to average over all possible pairs of 
particles 

n mixed (r \ _ J_ ^r(ri, ...,r i + r, ...,r N ) 

l N f N ) 

The asymptotic limit of the OBDM gives the Bose-Einstein condensate density 

= f ( 2 - 88 ) 
and the condensate fraction is obtained by calculating the asymptotic ratio 

oi (r) Nn 

lim = 2.89 

2.5.4 Two-body density matrix 

The two-body density matrix (TBDM) (^(^^(r^^^^^i)), depends on 4 vector 
arguments corresponding to destroying two particles at positions ri and r 2 and inserting 
them at positions r[ and r' 2 . The diagonal element ri = and r 2 = r 2 (see Eq. (1.8)) 
is called the pair distribution function. Generally, it depends on two vector arguments ri 
and r 2 . In a translationally invariant system (e.g. in a homogeneous gas) it is a function 
of (ri — r 2 ), that is a function of the distance between a pair of particles and two angles. 
We consider spherically averaged pair distribution function 

Let us explain now how this formula is implemented in Monte Carlo calculation. We 
make summation over all pairs of particles: 

.. N(N -l)JJ( ri -r a -r M R)|'dR 
92{) ~ n 2 L* /|^(R)| 2 dR nN /|^(R)| 2 dR 1 ' 

If we do a discretization of the coordinate space with spacing h and introduce function 
$h{z) which is one if z < h and zero otherwise (the distribution is obviously symmetric) 
one obtains the following expressions: 
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• In one-dimensional system: 

^) = (^E^(l^-^l)l) (2-92) 

In an uncorrelated system t^d^l) = 2h/L is constant and g2{z) = 1 — 1/N. The 
form of the pair distribution function depends on a dimensionality of the problem. 

• In a three-dimensional system the expression is 

Notice that the distance r enters explicitly in the expression of the pair distribution 
function, leading to larger numerical variance at small distances. 

2.5.5 Pure estimators and extrapolation technique 

In a VMC calculation one obtains a variational esimator for a given quantity (let it be 
denoted by an operator A), which corresponds to an average over the trial wave function 
i\) T : 

(4- = (2-94) 

Instead, the DMC method asymptotically provides a more precise mixed estimator 
given by 

(4- = (2.95) 

mm) 

Nonetheless, this type of average can differ from the inbiased ("pure") ground state 
average, which corresponds to the true quantum-mechanical equilibrium value at zero 
temperature 

r (0o|j|0o) 

{A puro = 2.96 
(0o|0o) 

The DMC method gives an exact result for the energy, as the mixed average of the 
local energy E\ oc = ip^HtpT coincides with the pure estimator. This can be easily seen by 
noticing that when (<f> acts on H, it gives exactly the ground state energy. 

From now on, we will demonstrate that averages of local diagonal operators can be 
calculated in a "pure" way. This means that the pair distribution function, potential 
energy and static structure factor can be estimated exactly. Local quantities are diagonal 
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in coordinate space (R|A|R/) = A(R)(R|R'). The "pure" average can be related to the 
mixed one in the following way 

. _ W^(R)j^l#) _ M(R)P(R)U ,,„- 



where P(R) is defined as 



P(R) = T^^ol^M • (2.98) 
^t(R) 



The asymptotic number of offsprings of a walker can give 0o(R)/Vt(R) [LKC74]. By 
"tagging" walkers one can identify, at any time, what is the parent walker. This forward 
walking method [BRL91, RBHL86, RR90, Run92] permits to sample pure averages. 

A more simple algorithm was devised by J. Boronat and J. Casulleras [CB95a], in which 
the explicit history should not be recorded and one operates with the actual values of an 
observable. This method is used in our calculations. 

Eq. (2.98) gives the number of descendants of a walker R for large times r — > oo. 
Practically it is enough to wait a sufficiently large, but a finite time T. One makes mea- 
surements of a local quantity for all of the walkers, but calculates the average after the 
time T, so that each walker was replicated according to the weight -P(R). 

An important example of a non-local quantity, for which no "pure" is known, is one- 
body density matrix (see Eq. (2.85)); this quantity deserves a special attention. We will 
explain an extrapolation technique, which can be applied for finding averages of non-local 
operators. It is also worth noticing that extrapolation can be used in estimating diagonal 
quantities, for example, pair correlation function. 

Adopting the notation 5ip for the difference between the trial wave function and ground- 
state wave function 

<t>o = ih + S*i>, (2-99) 
the ground-state average can be written as 

(i) purc = (0 o |i|0o) = (ih\A*fr) + 2(0o|i|#) + (5i(j\A\5<iP) (2.100) 

If 5i/j is small the second order term (5^\A\5xjj) can be neglected. After substitution 
(4>o\A\5ip) = (i[)t\A\<Po) — (i/)t\A\i/)t} the extrapolation formula turns into 

{A)p U re ~ 2(A) m i x — (A) var . (2.101) 

It is possible to write another extrapolation formula of the same order of accuracy: 



(2.102) 



2.5. Estimators for physical quantities 



47 



An applicability condition for the extrapolation method is that the expressions (2.101) 
and (2.102) yield the same final result. Nevertheless, in certain situations the use of the 
second formula can be advantageous, as it does not change the sign of the function, which 
is not always true in the case of the first formula, when the function is very close to 0. This 
can be useful in an estimation of essentially positive non-local quantities, for instance, the 
one-body density matrix. 



Chapter 3 

Ewald method for polytropic 

potentials 



3.1 Introduction 

The behavior of many-body systems is often governed by the long-range Coulomb potential 
between charged particles. Numerical simulations of such systems are usually performed 
by considering a finite number of particles in a cell with periodic boundary conditions. The 
correct estimation of the potential energy in such systems requires of a summation over all 
images created by the periodic boundary conditions. For long-range interaction potentials 
such direct summation either converges slowly or it is conditionally convergent, making 
its evaluation computationally cumbersome. Instead, the performance of the calculation 
can be greatly improved by using Ewald summation methods [Ewa21]. In these methods, 
the slowly convergent tail of the sum in the potential energy is represented by a rapidly 
convergent sum in momentum space. The method is named after Paul Peter Ewald who in 
his pioneering work dated almost a century ago calculated the electrostatic energy in ionic 
crystals (a detailed derivation for the Ewald sums for the Coulomb potential can be found 
in the work of de Leeuw et al. [dLPS80]). An alternative approach to deal with long-range 
systems is proposed by Smith [Smi94]. In his method, the Hamiltonian and equations of 
motion are derived using constraints on the velocities of particles. Instead, in the following 
we will stick to a standard model for the Hamiltonian and will consider ways to improve 
the convergence in the potential energy. 

For a good performance in simulations of large iV-particle systems, a number of modi- 
fied summation methods has been developed. Historically, the first efforts to enhance the 
Ewald method consisted in looking for appropriate truncation schemes, but all of them 
were strongly dependent on the system properties, in particular on the system size. Tab- 
ulations of precalculated terms in both real space and momentum space sums [SD76], as 
well as polynomial approximations of the involved functions [dLB47, BST66, Han73], were 
also proposed to look for a balance between calculation time and truncation errors. Never- 
theless, these approximate methods suffer from error accumulation in simulations of large 
systems, and do not allow for reducing the overall 0(N 2 ) complexity of the original Ewald 
summation. The work of Perram et al. [PPdL88] was the first to give a way to optimize 
the splitting of the interparticle potential between the long-range and short-range parts to 
yield a total complexity of 0(N 3 / 2 ). A special modification of the Ewald method called 
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Wolf summation [WKPE99, EGA10], based on a damping of the Fourier-transformed part 
of the sum, was posteriorly developed in order to render the original Ewald summation 
more efficient for non-periodic systems and large model sizes. The Ewald technique was 
also applied to develop the method of evaluation of electrostatic potential near the surfaces 
of ionic crystals [Par 75, Par 76]. 

Another way for improving the Ewald method is to perform fast Fourier transform 
(FFT) of a reciprocal space sum on a mesh. The oldest algorithm of this kind is the so- 
called Particle-Particle Particle-Mesh (P 3 M) method, invented and improved to the com- 
plexity O(N) in calculation of forces by Hockney and coauthors in the 70's [HE88, HGE73]. 
The P 3 M technique is based on a distribution of the charge density on a grid using a cer- 
tain smooth assignment function and then the discrete Poisson equation is solved using 
FFT. This algorithm appeared to be less complex to yield 0(N In N) with an appropriate 
choice of the free parameters. The P 3 M algorithm was recently improved by Ballenegger et 
al. [BCLH08] for calculation of energies, bringing, as claimed, the maximal precision in the 
energy by an optimization of the "influence" function (a substitution of the potential in the 
Fourier-transformed Poisson's equation). For a comprehensive introduction to Ewald- and 
mesh-based techniques we recommend to refer to the cited work of Ballenegger and coau- 
thors where special attention is paid to the estimation of both sum truncation-imposed and 
grid-imposed errors. The extension of this method, called Particle Mesh Ewald [DYP93] 
(PME), makes use of the analytical form of the sum in the reciprocal space and evalu- 
ates potentials via FFT instead of interpolating them as P 3 M does. Although PME is 
slightly more complex than the P 3 M algorithm, it is still 0(N In N) and allows to reduce 
significantly the memory expenses. Later Particle Mesh Ewald method was reformulated 
by Essmann et al. [EPB + 95], making use of cardinal i?-splines to interpolate structure 
factors. This approach, called Smooth Particle-Mesh Ewald (SPME) substantially im- 
proved the accuracy of PME with a comparable computational cost, as it still scales as 
0(N In N). SPME is also claimed to be applicable to potentials of the polytropic form 
l/|r| fc . In general, the conventional FFT-based approaches suffer from the severe fall- 
back of requiring equidistant particle positions. The invention of the variant of Fourier 
transform for nonequispaced nodes (NFFT) opened a path to overcome this shortcoming, 
while keeping the introduced errors below the specified target levels. The nonequispaced 
fast Fourier transform is currently considered as a promising means to improve the Ewald 
summation performance, with open code implementations available [KKP]. The early vari- 
ants of the NFFT algorithms are reviewed in the work of A. F. Ware [War98]; a general 
approach to the fast summation methods based on NFFT can be found in the article of 
G. Steidl [Ste98]. 

The most recent family of algorithms based on the Ewald approach are the tree-based 
algorithms, with the fast multipole method (FMM) being the most known and widely used 
among them. The algorithm, developed primarily by L. Greengard and V. Rokhlin [GR87], 
is based on the idea of keeping the direct summation of potentials or forces for the nearby 
atoms and approximating the interactions of the distant atoms by their multipole ex- 
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pansions. FMM offers the asymptotically fastest performance among the Ewald-related 
algorithms, being linear in N in most cases and not worse than O(NlnN) with explic- 
itly controlled accuracy. The FMM technique is naturally applicable to inhomogeneous 
and non-periodic systems, being also easy to parallelise since it is an entirely real-space 
summation. Since then the algorithm was significantly improved in efficiency, mostly by 
introducing new diagonal forms of translation operators [GR97]. However, FMM has an 
intrinsic shortcoming, when applied to molecular dynamics calculations, as the energy con- 
servation it brings is poor; the method per se is also rather cumbersome in implementation. 
Another group of methods, based on the multigrid methods of solving elliptic (in this par- 
ticular case - Poisson's) equations [BLdL98], was developed a decade ago [SD01]. These 
methods allow to preserve the scaling 0(N) and parallelization advantages of tree-based 
methods, as well as the applicability in simulations without PBC, being on the other hand 
satisfactorily energy-conserving and additionally accelerated on all length scales. An effi- 
cient realization of the multigrid method and its analysis may be found in the work of Sagui 
et al. [SPD04]. An advanced mesh Ewald technique, claimed to reduce significantly the 
computational costs of charge spreading in multigrid-based methods, was recently proposed 
by Y. Shan and coworkers [SKE+05]. 

A detailed comparison of the optimized 0(N 3 ^ 2 ) pure Ewald technique, FFT-based 
summations, and multipole-based methods was made by H. G. Petersen [Pet95] for systems 
with approximately uniform charge distributions, taking into account a possible parallel 
implementation. According to Petersen, the method of choice with a number of particles 
below 10 4 is the conventional Ewald summation, PME is preferable in the range N ~ 
10 4 — 10 5 , and the fast multipole method should overperform them with N > 10 5 . A 
more recent and ample review of FMM, P 3 M and pure Ewald methods by Pollock and 
Glosli [PG96], based partially on their own calculations, implies that P 3 M is faster than 
the Ewald summation already for 500 particles, although it is stressed that the other 
factors as the ease of the coding, the system geometry, as well as the code optimization can 
change the choice. We would also suggest a thorough survey of different Ewald summation 
techniques given in the work of Toukmaji and Board [TJ96]. 

An approach, alternative to using cubic periodic boundary conditions in a calculation 
of long-range interactions, called Isotropic Periodic Sum (IPS), was recently proposed by 
Wu and Brooks [WB05]. The main goal of their approach is to deal with long-range 
interactions, avoiding artificial correlations and anisotropy bias induced by a PBC-based 
summation in a cubic box. In this technique, only the interactions of a particle A with 
the others within a certain radius r c are taken into account (as in a plain cut-off scheme), 
and this spherical simulation zone is repeated in an infinite number of shifts by vectors 
r s h, such that |r s h| = 2N\r c \. Therefore, the particle A interacts not only with B (within 
the sphere radius), but also with all the images of B, occupying homogeneously the shells 
of radii 2A|r c |, centered in B. The subsequent integration and summation over the shells 
allows to obtain explicit expressions of forces and energies for a number of interactions of 
most physical interest, like electrostatic, Lennard- Jones and exponential potentials. The 
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method is known to yield a performance close to the one shown by the Ewald summation, 
but without imposing unwanted symmetry effects. 

Since its proposal, the Ewald method has been applied to a large number of physical 
problems, although mostly to systems with the Coulomb l/\r\ interaction potential. In a 
recent work by R. E. Johnson and S. Ranganathan [JR07], a generalized approach to Ewald 
summation is stated to obtain potential energy and forces for systems with a power-law, 
Yukawa potential and electronic bilayer systems. The Ewald method for two-dimensional 
systems with electrostatic interactions was developed by Parry [Par 75], but his technique 
appeared to be computationally inefficient. Spohr et al. [Spo97] studied a slab geometry 
by treating the simulation cell as a fully three-dimensional one with the conventional Ewald 
summation. Later on, a significant advance was made by Yeh and Berkowitz [YB99], as the 
authors managed to obtain the explicit correction term for the rigorous three-dimensional 
Ewald summation, that brings the results for a slab system in a satisfactory agreement 
with the 2D summation. The 2D Ewald technique was also applied by Wen Yang et 
al. [YJ06] to calculate the energy of Coulomb particles in a slab system with a uniformly 
charged surface. One of the first two-dimensional variants of the Ewald summation was 
presented in [GZW97], applied to the quasi-2D Stockmayer model with the potential l/|r| 3 . 
Recent applications to dipolar bosons in a 2D geometry have been made by C. Mora et 
al. [MP WO 7] and Xin Lu et al. [MPLW08]. On the other hand, the explicit forms of 
the Ewald sums for Yukawa interactions have been also reported: in 3D geometry, with 
partial periodic boundary conditions [SCOO, Maz07a], and in 2D geometry [Maz07b]. The 
general approach to the Ewald summation in quasi-two-dimensional systems with power- 
law potentials and results for several values of power factor are given in [MazlO]. The 
Ewald method can also be useful even applied to fast decaying power-law potentials. For 
instance, the Ewald formalism was developed in [KI89] for the dispersion interaction l/|r| 6 
and later for the Lennard- Jones potential by W.-Z. Ou-Yang et al. [OYLS + 05]. Also, Shirts 
et al. in their recent work [SMCP07] argue the need for taking into account the effects of 
cutoffs in molecular dispersion interactions due to a Lennard- Jones potential, especially in 
non-isotropic and inhomogeneous media. The authors developed two formalisms for the 
estimation of these cutoff errors in binding free energy of macromolecular systems, which 
can in principle be extended to the other observables. However, it is claimed that the 
adequate implementation of the Ewald summation for this kind of systems may render 
their corrections unnecessary by mostly eliminating the cutoff-dependent behavior. 

In this Chapter, we report explicit expressions of the Ewald sums for the general case 
of particles interacting via a l/|r| fc polytropic potential and in 3D, 2D, and ID geome- 
tries [OAB12]. The closed derivation of these sums is given, with special attention being 
paid to conditionally convergent potentials. One of the difficulties of the derivation is that 
different terms have to be considered in the cases of short-range, long-range or "marginal" 
potentials. In the case of a short-range interaction, the original slowly convergent sum is 
represented as a linear combination of two rapidly convergent ones. For a long-range in- 
teraction, the condition of charge neutrality in the simulation cell is shown to be necessary 



3.2. Ewald method in 3D geometry 



53 



to make the energy absolutely convergent within the considered scheme. The introduction 
of a uniform neutralizing charged background (jellium) , as a particular case of a charge- 
neutral system, is also discussed. The explicit forms of the Ewald sums are reported for 
a jellium system and for an arbitrary polytropic potential. We explicitly calculate the 
expressions for physically relevant interactions as Coulomb, dipole-dipole, and Lennard- 
Jones potentials. Finally, we have extended the Ewald sums to the case of a noncubic 
simulation cell, that could be useful in simulations of hexagonal closed packed (hep) and 
two-dimensional triangular solids. In addition, the general derivation path given in this 
work may be used to obtain the forms of Ewald sums for other interaction potentials. 

The computational efficiency is another important issue of the practical implementation 
of the method. In fact, one needs to choose correctly a free parameter, appearing in the 
integral representation of the sums, and to decide which number of terms should be kept in 
spatial and momentum sums in order to reach the required accuracy. The choice of these 
three parameters affects the difference between the calculated result and the exact one as 
well as the calculation complexity. Therefore, a certain optimization of the parameters 
is always required. In this Chapter, this optimization process is formalized and it is 
shown that following the described procedure the overall computation time is significantly 
reduced. The accuracy of the result is shown to be kept under control, with the only cost 
of a preliminary benchmark calculation. 

The rest of the chapter is organized as follows. In Section 3.2, we formulate the problem, 
develop the general Ewald approach and report explicit expressions for the Ewald sums for 
a polytropic potential in a three-dimensional cubic simulation cell. Sections 3.3 and 3.4 
contain derivations of the Ewald sums in two-dimensional and one-dimensional geometries, 
respectively. In Section 3.5, the case of a simulation cell with different side lengths is 
considered for three- and two-dimensional systems. The final general expressions and their 
particularization to the most physically relevant cases are presented in Section 3.7. The 
practical algorithm for the parameter optimization and an actual application of the Ewald 
method is discussed in Section 3.8. Summary and conclusions are drawn in Section 3.9. 

3.2 Ewald method in 3D geometry 
3.2.1 Basic assumptions and initial sums 

We consider a system of N particles inside a cubic simulation cell of size L with periodic 
boundary conditions. Thus, each particle with coordinates r in the initial cell has an 
infinite number of images r + nL in the adjacent cells. The total potential energy is 
estimated by 




N N 




nez 3 



i=\ 3=1 
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where <p(r) is the interparticle potential, = Vi — rj, and the prime in the first sum means 
that the summation over an integer vector n must be done omitting the term n = when 
i =j. 

3.2.2 Analytic derivation 

In many physical situations, the interaction potential between two particles i and j has 
the power-law form qiqj/\r\ k with positive k and g,, qj being the generalized charges of the 
particles. This sort of interaction is generally referred to as polytropic potential. 

First, let us consider the case of short-range potentials, k < 3. As we will see later, 
the potentials corresponding to k > 3 give a similar result. For k < 3, the right-hand 
part of Eq. (3.1) diverges and it can be made convergent only if the restriction of charge 
neutrality is required, i.e., when X^=i q% — 0- It has also been shown [FF96] that for a pure 
electrostatic interaction the total energy (3.1) can be conditionally convergent even in a 
neutral simulation cell because of a higher multipole contribution. The energy and forces 
are therefore dependent on the order of summation, which can also be implicitly set by 
a choice of a convergence factor. The ambiguity usually appears in a form of a constant 
or a position-dependent term, vanishing in the limit L — > oo. Hence, the preference 
in one or another factor should be dictated either by physical properties of a particular 
system or by arguments regarding rates of convergence to the thermodynamic limit. For a 
general discussion on the convergence issues appearing in periodic boundary conditions, see 
Ref. [MP95]. The main idea of the Ewald summation technique in the approach proposed 
by de Leeuw, Perram, and Smith [dLPS80] is to multiply each component of the sum by 

2 

the dimensionless factor e~ sn , with s > being a dimensionless regularizing parameter, 
making the sum absolutely convergent. Then, the limit s — > is taken, so that the 
singularity in the initial sum (3.1) can be explicitly separated into a term depending only 
on s, that finally can be cancelled due to the charge neutrality condition. We take a similar 
multiplier c(n, r, s) = e^ s ^ n+r ^ yielding the same rate of convergence (since < r < 1 in 
units of L). As the sum, multiplied by c, is invariant to an arbitrary substitution r — > n+r, 
the chosen convergence factor allows to preserve the periodicity of the potential in order 
to avoid any possible artefacts in the final results. 

For the sake of clearness of the derivation, it is convenient to use reduced length units, 
that is to use the size of the box L as unity of length and substitute by r^L. From now 
on, and for simplicity, we use the notation r for r»j and, in case of possible ambiguity, we 
will stick to the standard notation r^. Also, we rewrite the potential energy by splitting 
the total sum (3.1) into two terms: ioi (the sum of the interactions between a particle with 
all the other particles in the box), and Too (the sum of the interaction of a particle with 
its own images, comprised of the components i = j in Eq. (3.1)). Explicitly, 




(3.2) 
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with 



'01 



E 



E 



i<i<j<N \ r V + n l 



1 e" sn2 



(3.3) 



J oo = 2 E V^ g - (3 - 4) 

z nSZ3\0 U t=l 

where the shorthand notation n = |n| is used. 

First, let us focus on the /oi term, which we rewrite as 

J oi = E QiQji J ( r ^ s ) > ( 3 - 5 ) 

l<i<j<N 

where we have defined the "screened" interaction potential ip(r, s) = J2 n e _ ' s ' r+n ' 2 /|r + n| fc , 
extended from a single cell to the whole coordinate space. Since the total potential energy 
consists of a sum of pair interaction components, we may consider a single pair without 
any loss of generality. 

Let us apply the equation 

roc „ 

x~ 2s = -f- / t s - l e~ tx dt , (3.6) 

r(s) Jo 

representing the definition of the gamma-function, to the polytropic potential \r + n\~ k . 
Then the function ip may be represented in an integral form, 

tf(r, a) = jH ft" 1 £ e-^V^I 2 dt . (3.7) 

We expect that the integral (3.7) contains a singularity that will be located in the 
vicinity of zero. Therefore, we split this integral into two domains [0, a 2 ] and [a 2 , oo), the 
corresponding integrals being denoted as i/j^ and il>mf, where a is some arbitrary positive 
constant, 

ij;(r, s) = ^ fin (r, s) + ip in{ (r, s) . (3.8) 
In the following, we analyze the two terms of the previous sum (3.8). 
1. The explicit analytical form of the term Vw( r ; s) can be found 



Mr, s) = -^-E r ti-*e-W-W dt = £ 

1 {2) n Ja n 



e-*l r+n l a r(|,a 2 


r + n 


2 ) 


\r + n 


k y 


(!) 





(3.9) 

where T(a, z) is the incomplete gamma function. From the large distance asymptotic 
expansion of this function, one obtains that the above lattice sum is absolutely and 
uniformly convergent if s > and a > 0. Therefore, one may simply take the limit 
of vanishing screening s — > 0, 

^0 1 rj^+np) 

iMr, ») -* E | r + n[fc • ( 3 - 10 ) 
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2. The calculation of ips n (r, s) is done by making a separate analysis of the n = case, 



^(r )S ) = C/°(M) + C =0 (v) 



Explicitly, 



fin I 



3 

7T2 



r,s 



E 



3 

7T2 



(t + s) 
df 



r exp 

2 



-7r 2 m 2 



t + s 



+ iTx'vmr 



df 



a 2 ^-1 



r(|)^o + 

where we have used the Jacobi transformation [Jac21, WW] 

'7T\ 3 / 2 



(3.11) 

(3.12) 
(3.13) 



applied to 



-s|ra+r| 



exp[— 7r 2 m 2 / s + 2iiimr] for n, m e Z 3 



exp[— s|n + r| 2 — t|n + r| 2 ] = exp[— (s + t)\n + r\ 2 }. 



(3.14) 



(3.15) 



We evaluate the integral ip^f°(r, s) by the following analysis. Consider separately 
the following factor of the integrated expression from (3.12) 



M 



exp 



t+s 



(t + S 



(3.16) 



It is clearly continuous and bounded on (0, +oo) as a function of (t + s), also notice 
that t fc / 2_1 is absolutely integrable on (0, a 2 ) for k > 0. In accordance with the 
standard convergence test for improper integrals, the integral ip^f°(r,s) converges 
absolutely and uniformly with s being considered as a parameter. Then, the limit 
s — > may be carried out and the integral becomes 



3 

7T2 



r,s\ 



E 



^2irimr 



r(k) 



fc-5 

t 2 exp 



E 



7T2 cos(27rmr) k _ 3 



n 2 m 2 



' -K 2 m 2 ^ 



df 



r(|) 



a" 



(3.17) 
(3.18) 



The function E n (z) is the exponential integral function, and we have cancelled the 
imaginary part of the sum (3.17) by grouping the pairs with n and —n. 

Now, we analyze the second term of ipan{i", s), 



m=0 1 

fin 1 
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7T2 



r,s 



« 2 ft" 1 



r(f)^o (t + s) 



dt . 



(3.19) 
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In terms of a new variable v = s/(t + s), 

CT>, s) = ^rf] (1 -^~V 3)/2 dt, . (3.20) 

The integration of ^^ =0 (r, s) for a l/|r| fc interaction has to be carefully analyzed as 
a function of k: 1 < k < 3, long-range potential; = 3, marginal case; and k > 3, 
short-range potential. 

(a) Suppose 1 < fc < 3. The resulting integral, 

CT°(r. s) = CTV, s ) = ^rJ 1 /l2 , (1 ~^"^ (fc - 3)/2 ^ (3.21) 

1 1.2/ Js /( a + s ) f 2 

may be given explicitly in terms of incomplete beta- and incomplete gamma- 
functions. Expanding the resulting function for small s, 

It is easily seen, that the only divergent term in the expansion (3.22) is the first 
one, which we define as 

(3.23) 

We remind that the choice of a convergence factor (that explicitly affects the 
summation order) may in principle lead to additional contributions in the total 
energy if the convergence of the sum is conditional (like for a charge-neutral 
cell of Coulomb particles with non-zero total dipole moment). In the origi- 
nal derivation of de Leeuw et al. [dLPS80], the factor exp(— sn 2 ) results in an 
additional dipole-like component in ip^ , which breaks the periodicity of the 
potential and therefore complicates its use in simulations with periodic bound- 
ary conditions. Moreover, this procedure [dLPS80] yields a nonvanishing dipole 
term exclusively for k = 1 in 3D geometry, with the rest of the sums remain- 
ing unchanged. From our point of view, this discontinuity points out to an 
nonphysical character of the dipole term appearing in the case of the Coulomb 
potential. Nevertheless, in a number of studies [FF96, MP95] it is considered 
as a first order correction when the convergence to the thermodynamic limit 
is analyzed. The mere fact that the results for the two different convergence 
multipliers coincide when k > 1 is a consequence of the absolute convergence of 
the higher multipole contributions in this case. 



3- k 



+ 



27T2a 



fc-3 



(fc-3)r 



O(s) 



(3.22) 



S(s) 



fc-3 „ 

S — 2vrr 



3 - k 
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(b) Suppose k = 3. In this marginal case, the expression (3.21) may be integrated 
directly to yield the following logarithmic dependence 



m=0 , 
fin 1 



r,s 



3 

7T-" 



-2as - 2a 3 



[a* + sj2 
that close to s = expands as 



+ ln(s + 2a 2 + 2aVs + a 2 ) - In s I (3.24) 



m=0 1 
fin 1 



r, s) = -2ir In s - 47r + 47r ln(2a) + 0(s In s) 



(3.25) 



with the diverging term 



3(s) 



-27rln s . 



(3.26) 



k 

(c) Consider the remaining option k > 3. In this case, (1 — v)2 is bounded from 
above and (k — l)/2 > 1. It means that the integral converges absolutely and 
the only finite contribution to the integral comes from the first (constant) term 
of the integral expansion for small s, 



tf, 



m=0 , 

fin I 



27T2 a 



k-3 



r, s 



[k-3)T 



(3.27) 



The second term of the total potential energy, loo (3-2) can be derived in a similar form 
to the first one. The procedure to find the form of ip(r, s) is repeated here with = 0, 
hence the results are obtained straightforwardly via (3.10), (3.18), (3.25) and (3.27), 



'00 



E% 2 

i=l 



r 



1 r(|,a 2 n 2 ) 

rk\ 2.^ „k + 2^ 



3 

7T-" 



2> n 



.a k ~ 3 E k -i 
r k) — 



Q- 



r(f + i) 



m=0 , 
fin 1 



r, s 



(3.28) 



with the term ip^ °(r, s) depending on the potential parameter k via (3.22), (3.25) or 
(3.27). 

Putting all together, the potential energy can be written in a more compact form as, 



U 



^jfcC^Ol + ^oo ) 



1 1^1 1 

Tk E ngrtto/L) + — £ + Tk E mS(s) + 7771 E rf^M , (3-29) 



L k 



2L k 



i=l 



L k 



i<j 



2L k 



with the generalized potential, 



ifj(r) = R(n, r) + K(m, r) + d . 

n m^0 



(3.30) 
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A constant shift in the definition of ifj is introduced to satisfy by the property / cell ipdr = 0, 
convenient for a proper treatment of the background contributions (see Appendix A). The 
functions entering in Eq. (3.30) are defined as 



r(|,« 2 


r + n 


2 ) 


r(|)|r + n 


k 



R(n,r) = v jy ' ■ ' (3.31) 



K(m, r) = x(m) cos(27rmr) , (3.32) 

(3.33) 



with 



3 



fe-3 / J2^2 S 



2 

The explicit form of the function S(s) depends on the k value 



fe-3 „ 

s— 2vrr 



3-fc 



if < 3 (singular term) 
S{ s ) — { -27rlns if k = 3 (singular term) , (3.35) 

if k > 3 



and the term £ depends only on the choice of a, 



£ = £ p(n) + X *M + Ci + C 2 , (3.36) 



with 



and x(n) defined in Eq. (3.34). The constants C\ and G% are explicitly, 



27r^q fc - 3 

TFT 



if fc^3 



d = <( (*-3)r[$j ' (3.38) 

-4vr + 4vr ln(2a) if fc = 3 

C2 = 7^— (3.39) 

r(f + i) 1 ; 

(3.40) 

3.2.3 Removing singularities for k < 3 

The diverging part t/ s (containing a singularity) of the total potential energy equals to 

Us = ^ E fl^to + £ <7^( s ) = 2T f E ^ ( 3 - 41 ) 

i<j % \ i / 

and vanishes, if the charge neutrality condition Yli Qi — is taken. 
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Consider now a charge-neutral system with a neutralizing background consisting of a 
large number of identical uniformly distributed particles of the opposite charge (the "jel- 
lium" model). We denote the numbers of negatively charged particles g_ and positively 
charged (background) particles q + as N_ and N + , respectively. By imposing charge neu- 
trality, q + = — [i\C/iV + ]g„, with N the total number of particles, N = iV_ + N + . 

The potential energy for the jellium model can be written as 

U =JkJ2 wMra/L) + N - qi + L ? +q2+ Z ■ (3-42) 

i<j 

The second term in Eq. (3.42) has a component proportional to N + q+. Note that the 
negative charges g_ and their number N_ is defined by the problem and therefore fixed. 
Hence, in the limit N + — > oo, this term cancels N + q^_ = {N"^q 2 _)/N + — > 0, and therefore 
this background contribution may be eliminated to yield 

^^=^- (3-3) 
Concerning the first term of Eq. (3.42), let us split it into three pieces, 

7 E MMr) = j (S- + 2S-+ + S ++ ) , (3.44) 

l<i<j<N L 

where the first sum corresponds to the interaction between the negative charges 

S—= E QiQMr), (3-45) 

l<i<j<JV_ 

the second sum is the interaction of the negatively charged particles with the positive 
charges of the background 

N— N + +N_ 

5- + = E E QiQMr), (3.46) 

4 = 1 j = l + N- 

and the third one is the interaction between the background charges 

s ++ = m^ir) . (3.47) 

l+N-<i<j<N++N- 

The last two terms S ^ and S ++ are easily shown to be zero in the limit N + — > oo as a 

consequence of the zero value of the integral of ip over the simulation cell (see Appendix 
A). 

With the above considerations we can finally write the expression for the potential 
energy within the jellium model as 

U icl = %T.^/L) + ^ (3.48) 

i<j 
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In the more general case of different charges in a charge-neutral simulation cell (with 
a long-range potential) or a system with an arbitrary short-range potential the potential 
energy is given by 

1 T N a 2 

U scn = JT E wMnj/L) + • (3-49) 

i<j 

A certain analytical conversion of the sum in the reciprocal space is also possible in 
order to sum it up faster. Expanding the sum that defines K(m, r) (3.32), one can simplify 
it in the following way, 

E<Mj E K(m,r) = x{m)^q i q j cos{2'Kmr) 

i<j mj^O m^O i<j 

= 2 E x ( m ) E QiQj [cos(27rmrj) cos(27rmr,,) + sin(27rmr,;) sin(2irmrj)] 

m^O i,j 

-^E% 2 E x M 



- E 



qj exp(2irimrj] 



i rra^O 



In this form, the sum over all pairs of particles in the reciprocal space is represented as 
a single sum over particles and thus it scales as O(N) instead of 0(N 2 ). Notice that the 
number of prefactors x(m) and exponents in the sum depends on a chosen cutoff, which in 
general also might depend on N, making the overall complexity of the fc-space grow. Naive 
schemes with a and the cutoff not depending on N do not take into account the interplay 
between the r-space and fc-space sum complexities, thus leaving at least 0(N 2 ) in one 
of them. Nevertheless, as we show later, optimization with a and cutoff depending on N 
gives a best total complexity of 0(N 3 ^ 2 ). An alternative method to sum up the momentum 
space part is to use Fast Fourier transform-based techniques (like PME), which is fast as 
0(iVln N). 

The last term in Eq. (3.50) cancels the x(m) component of £. Introduce the notation, 

$(r) =Y,R(n,r) + C 1 (3.51) 

n 

| = E P( n ) + Ci + C 2 (3.52) 

5equai(m) = <?_ E exp(2nimr j /L) (3.53) 

j 

S q (m) = qj exp(27rimr,/L) , (3.54) 

3 

where S^quai is used when the system of equally charged particles g_ is considered. Within 
this notation the potential energy may be rewritten in the following forms, which are more 
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efficient for numerical implementation, 

^ jel = f| E#W^) + oT* £ M™)|£ cqU ai(m)| 2 + ^1 (3.55) 

f/ gcn = TkYtHlM^/L) + — £ x(m)|^(m)| 2 + , (3-56) 
with rj, rjj in the original length units. 



3.2.4 Short-range potentials and the marginal case 

In case of a short-range interaction (k > 3), the potential energy does not diverge, which is 
clear from the form of the singular term S'(s)(see Eq. 3.35). Hence, there is no need to add 
a neutralizing background and, even more, the background must be necessarily excluded 
since it leads to a divergence in the energy. This is easily seen by considering the potential 
energy of the background separately 



/■cell d r 

U bg = cj — , (3.57) 



that contains a singularity in zero. The expression for the potential energy is simply equal 
to Eq. (3.49), 

1 T N a 2 

i<j 

When k — 3 (marginal case), both ultraviolet (coming from short-range contributions) 
and infrared (coming from long-range contributions) divergences arise in zero for the back- 
ground as well as in the vicinity of infinity (the logarithmic divergence in the energy of 
negative charges). The only coherent model here is a plain "quasi-neutral" gas consisting 
of a mixture of a finite number of charges per box with the constraint X) = 0, i.e., with 
the positive background excluded. 



3.3 Ewald method in 2D geometry 
3.3.1 General notes for lower dimensions 

The Ewald sums can be extended to two-dimensional (2D) systems interacting through 
polytropic potentials. The difference with the 3D case comes from a different form of 
the Jacobi imaginary transformation for the Jacobi ^-functions [its 3D form is given in 
Eq. (3.14)]. 

The "third" Jacobi ^-function 8 3 (z,t) is defined as 

+oo 

3 (z| r ) = Yl e i7TTm2 e 2miz , (3.59) 

m=—oo 
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and satisfies the Jacobi imaginary transformation, 

9 3 (z\r) = (-ir)- 1 /^' 2 ^ e 3 (zr'\r') , (3.60) 

with t' = —1/t. Under the change of variables, z = irr and r = ivr/s, the ^-function 
becomes a Gaussian, which is the relevant function for performing the Ewald sums, 

J- e -<r+n)* = (^5)1/2 e -" 2m2/s e 27rimr . (3.61) 

n=— oo m=— oo 

This expression will be used later, in the derivation of the Ewald sum in one-dimensional 
systems. Equation (3.61) may be easily generalized to the 2D geometry, 

]T e" s|r+n|2 = (tt/s) e-" 2m2/ V™ . (3.62) 

n m 

Comparing this result for 2D with its ID (3.61) and 3D(3.14) counterparts one finds that 
the dimensionality D affects only the constant multiplier as (tt/s) d ^ 2 . 



3.3.2 Analytic derivation 

The analytical derivation of the Ewald sum in 2D proceeds similarly to the one already 
presented for 3D. Equations from (3.2) to (3.11) are also valid here because their derivation 
is done without explicit reference to the dimensionality of the problem. In particular, the 
integral -?/>mf (r, s) converges absolutely and to the same value 

1 ^ r(|,a 2 |r + n| 2 ) 

ipmf(r,s) — 4 ,, s > — —. r . . (3.63) 

r(|) n \r + n\ k v ; 

We make the same decomposition of the integral ipa n (r, s) as in 3D, 

^(r, s) = ^Sf»(r, s) + ^=°(r, s) , (3.64) 

with 

where the two-dimensional variant of the Jacobi transformation (3.62) is used. The differ- 
ence between the pair of equations (3.65, 3.66) and their three-dimensional analogues (3.12, 
3.13) relies in a substitution of the 3D factor (%/ (t + s)) 3 / 2 by the 2D one ir/ (t + s). 
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(3.65) 
(3.66) 
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First, we consider the term ip™f°{r, s). Following the same analysis as for its 3D 
counterpart, it can be shown that this parametric integral also converges absolutely. It 
yields 



7T 
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2mmr 



t 2 exp 



7r 2 m 2 



dt 



\- 7rcos(27rmr) fc _ 2 

, r -) 2 



ix 2 m 2 ^ 



a* 



(3.67) 



The modification of the integral ^fkT ^ s ^ ess straightforward, since it requires specific 
integrations and expansions in the series for small s. Namely, we have to evaluate the 
integral 



m=0 

fin 



7T 



(1 



V 2 



1)1 



(3.68) 



which is the 2D equivalent of Eq. (3.20). 

In the following, we consider separately the cases of long-range potential (1 < k < 2), 
marginal interaction (k = 2) and short-range potential (k > 2). 

1. 1 < k < 2. As in 3D, the integral can be found analytically via the incomplete beta- 
and incomplete gamma-function with known series expansions for small s. Omitting 
these unnecessary intermediate expressions, we give the final expansion for V ; ft n =0 ; 

^ 2 27ra k ~ 2 



1>: 
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fin 
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sin (^) r (|) (k - 2)r 
The first term of the expansion, 
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(3.69) 



(3.70) 



clearly diverges when s — > 0. Similarly to the 3D case, this term is cancelled in a 
charge- neutral cell and hence, 

2na k ~ 2 



ih m=0 
^fin 



(A;-2)r 



(3.71) 



2. k = 2. The integration of Eq. (3.68) is performed to yield in the limit s — > a 
marginal logarithmic dependence, 



m=0 

fin 



— 7r In s + 2ir In o + 0(s In s) . 



(3.72) 



As for the 3D geometry, the jellium model is inapplicable in this particular case 
since the energy of the continuous background diverges. Nonetheless the diverging 
component 

S(s) = -7rlns (3.73) 
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can be removed if we consider a charge-neutral system with a finite number of charges. 
In this case, 

^™=° = 2vrlna . (3.74) 

3. k > 2. The integral (3.68) can be evaluated by taking s = 0, since its convergence is 
absolute, 



(fc-2)r 



The second potential energy component, I 00 (3.2), is calculated as in the 3D case. The 
result for 2D is 
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2/ / m^O 



loo(s) = E% 2 (^T (0, S )+^nf(0, S )-Cf(0, S )) (3.76) 
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3.3.3 Final expressions 

With respect to the 3D case, the changes in the 2D Ewald sum appear in those terms in 
which the Jacobi transformation is used, that is in x(n) and C±, 

, s ira k ~ 2 _ ( ix 2 m 2 \ . „„. 

C, = Cr° (3.78) 

The other terms, namely R(r,n), p(n) and C2, are not affected by dimensionality and 
may be taken directly from the previous section. 

Within the jellium model for a long-range potential (k < 2), the Ewald sum is given by 

^ = |E^/ i ) + ^- ( 3 - 79 ) 

A more general form, applicable to any system with a short-range potential {k > 2), a 
charge- neutral system with long-range interaction (k < 2), or a marginal (A; = 2) potential 
is expressed as 

1 £ * 

^ gen = 71 £ qwttm/L) + £ ^ 2 • (3.80) 

i<j i=l 

In the same way as for the 3D systems we can modify the sum in the reciprocal space, 
and with the same notations (3.51) - (3.54) (p, R and the constants C±, C2 are the new 
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ones, corresponding to 2D case) the potential energy may be given by 

a 2 1 ~ Nn 2 ~ 

U icl = jtE^j/L) + — £ x(m)|S cqual | 2 + ^ (3.81) 

i<j m^O 

= -^Em^/L) + ^ E *(™)\S q \ 2 + §f f • (3.82) 

3.4 Ewald method in ID geometry 

As it has been commented before for the 2D case, the differences due to dimensionality 
are caused by the form of the Jacobi imaginary transformation. In the derivation for ID, 
one needs the following ones 
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x~ 2s = -f- / " t^V** 2 dt (3.83) 
Vis) Jo 



+oo 

J2 e~ sn2 = (vr/s) 1 / 2 J2 e~ n2m2/s (3.84) 

n=— oo m=—oo 
+oo +oo 



m=—oo 



Similarly to what discussed in the previous section, the only terms to be changed are those 
where the Jacobi transformation is used, namely ip™^° (in /oi in a radial- dependent form, 
in Iq Q for r = 0). The difference arises from a different power exponent (1/2) in (3.84) and 
(3.85), that is in (3.18) k has to be substituted by k + 2 (and 7r 3 / 2 - by vr 1 / 2 , respectively), 
yielding 

1/2 27rimr / 2 2\ 

As far as the term ^fin -0 is concerned, we should perform a simple integration and do 
a series expansion for small s, 



k-l 



The estimation of this integral depends on the k value. In the following, we detail this 
analysis. 

1. k — 1, the marginal case, 

Cr° = In 8 - 2 + 2 ln(2a)) + 0(s) . (3.88) 

As before, we keep only the constant term, considering the diverging term absent 
due to the charge neutrality condition. Therefore, with T(l/2) = ^Jlx one has 

CT° = -2 + 2hn» (3.89) 
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2. k > 1, the short-range potential, 



k 

2 



VC = ° = ■ I3T + + 0{s^ k - 1)/2 In s) . (3.90) 



In the limit s — >■ 0, it yields 

27r l/2 k-1 

resembling the 3D result (3.27), with the change k — > k + 2 (except in the T term) 
and 7r 3/2 7r 1/2 . 

The final result for the one-dimensional Ewald summation reads 

^(r) = J2 R ( n i r )+ E%> r )) + C i ( 3 - 92 ) 
£ = E Pfa) + E *H + Ci + C 2 , (3.93) 

where Ci = V'/in * s taken from the expressions (3.89) (if k — 1) or (3.91) (if A; > 1). 

For k — 1, the only consistent system is the charge-neutral one with a finite number of 
particles. In this case and for a short-range potential {k > 1) one the potential energy is 
given by 

1 T N a 2 

U scn = Z E QiQM r v/ L ) + ^fgr* • ( 3 - 94 ) 

Although the Ewald method is applicable to one-dimensional problems, there is a direct 
way to calculate the sums for polytropic potentials 

-i n=+oo -i 

U = 71 E 7" • (3- 95 ) 
L k n t^oo \r + n\ k 



For k > 1, this sum can be represented as a linear combination of the Hurwitz zeta 
functions, 

1 +°° i I 

71 E 7^^ = TsWr) + m ~ r)) ■ (3.96) 
L K n r7 OC) |r + n\ K L K 

In particular, for k = 2 the sum converts into a familiar expression used in the Calogero- 
Sutherland model [Sut71, AGLS06], 

1 +°° 1 2 

77 E .— = 2/ \ • (3-97) 

L 2 n £^ 00 \r + n\ 2 L 2 sin 2 (7rr) v ; 

Notice that the sum (3.96) may be expressed in terms of trigonometric functions only for 
even values of k via {k — 2) times differentiation of Eq. (3.97). Anyway, the possibility to 
find exact expressions for infinite sums in ID suggests that the use of the Ewald method 
might not be needed, but we keep it as a possibly useful mathematical relation and for 
completeness. 
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3.5 Generalizations to non-cubic simulation cells 



3.5.1 3D case 



A special and interesting situation arises if we consider a simulation cell in a more general 
way, as a rectangular box with different side lengths (L x , L y , L z in the corresponding di- 
mensions) . The need to deal with a box of unequal size lengths may occur in the simulation 
of a solid with a noncubic lattice (the simplest examples include a hexagonal closed packed 
crystal in 3D geometry), since the lattice vectors n in the sum over images on (3.1) are 
no longer orthogonal. Focusing our analysis to a 3D geometry, the potential energy is now 
given by 



U 



2 ^ 

n a ez 3 



N N 

EE^ 

i=\ j=i 



L n r ) 



(3.98) 



with n r = (n x L x + n y L y + n z L z )/L , n x>y>z being integer vectors along the corresponding 
axis x, y, z. We have introduced the geometric average Lq = {L x L y L z ) 1 ^ and we will 
use reduced L Q units for r^, and hence will be adimensional. Repeating the standard 
procedure, we multiply the potential energy by a Gaussian term exp(— s\n r + r\ 2 ) and, at 
the end, we take the limit s — > 0, separating the converging part, if present. We group 
separately the interaction with images of other particles /oi an d the interaction of a particle 
with its own images Iqq, 



1 



01 



( 00 J 



(3.99) 



where 
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'00 
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qiqje s \ n r+ r ijV 



1 s-^ e 

^ ngZ 3 \0 I 



l<i<j<N \' l J 
-s\n 2 r \ N 

— ft E <ii ■ 



i=l 



n. 



(3.100) 
(3.101) 



Comparing the relations (3.99) - (3.101) to the cubic case (3.2) - (3.4), one notices that 
these relations remain unchanged if n is formally substituted by n r , and the constant 
coefficient 1/L k is replaced by 1/Lq. Therefore, all the results found without the Jacobi 
transformation (3.14) remain the same with n r instead of n. In particular, Eq. (3.10) 
transforms into the following 
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(3.102) 
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The Jacobi transformation (3.14) in a noncubic box has the following form 

-s\n r +r\ 2 _ TT ^sjnjLj/ Lp+n) 2 

i=x,y,z Tn 

\ 1/21 / 



n. 



n 



TT 



i=x,y,z 



s(Li/L 



n E ex p 

i=x,y,z TUi 



n 2 m 2 



s(Li/L )\ 



/sf /2 exp(-7r 2 |m fc | 2 /s) exp(27rim fc r) 



exp(2'Kim i r i L /L i ) 
(3.103) 



with m,k = m x Lo/L x + m y Lo/L y + m z Lo/L z the normalized displacement vector in mo- 
mentum space. The last equation is obtained from the original expression (3.14) by a 
formal substitution of the vector m by m^. 

In order to calculate ^/v we first modify Eq. (3.15), 



exp[— s\n r + r\ 2 — t\n r 



exp 



t)\n r 



then insert it into the relation (3.103), and finally separate the summand n = 0, 
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o {t + sy. 
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exp 



-7r m\ 
t + s 



+ 2nim k r 



dt 
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7T2 



r(|)io (t + s )i 



(3.104) 



dt 



, fin ^« =U - (3-105) 

The subsequent derivation follows exactly the derivation for a cubic box, with the change 
of n by n r and m by for sums in the real and momentum spaces, respectively. The 
final result for a 3D system in a noncubic box can be summarized as follows 



u 



E 



T(k/2,a 2 \n r + r 



TT2a k 3 cos(27rm,fcr' 



r(fc/2)|ra r 



r(fc/2,a 2 |n r | 
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T(k/2) 
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7r|m fc | 
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n 5o T(fc/2)|n r 



fc-3 
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5o r (^/2) 



7r 2 1 | 



(3.106) 
(3.107) 

(3.108) 



with the constants C\ and C 2 defined in (3.38) and (3.39). As it was done in the cubic 
box, the potential energy may also be given with the momentum space sum (linear in N). 
Applying the definitions, similar to Eqs (3.51) - (3.54), 





= Y^R(n r ,r) + C 1 


(3.109) 




= PK) +C! + C 2 


(3.110) 








S e qual( rn k) 


= <?- E exp(2vrim fc r i /L) 

3 


(3.111) 


S q (m k ) 


= ^]g i exp(27rim fc r j /L) , 


(3.112) 
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the potential energy for a one-component jellium model converts into 

a 2 1 ~ Na 2 ~ 

U ]cl = E^fa/Lo) + 771 E x(m fc )|5 cqual (m fc )| 2 + ^ , (3.113) 

with a natural extension to the general case 

U" cn = iSft^Cry/Lo) + 7^ E x(m fe )|4(m fc )| 2 + |-f f • (3.114) 

Note that the formulas are derived for an orthogonal basis set. Still many triclinic lattices 
can be sampled in a similar form. In such cases the crystal is constructed not by translating 
the smallest-volume unit cell along non-orthogonal vectors but rather by translating a 
pseudo-unit cell of size L x x L y x L z containing more atoms along orthogonal directions. 
For example hep crystal can be summed in this way. Nonetheless, as the pseudo-elementary 
cell technique may be inconvenient in application to triclinic lattices, we would suggest the 
reader to rely on a reciprocal lattice technique (see [AM76]). 



3.5.2 2D case 

The generalization of the formulas found in a square 2D geometry to a rectangular simu- 
lation box comes in a similar manner. It is sufficient to take the resulting expressions for 
the two-dimensional problem (3.63) and (3.67), and to perform the necessary substitutions 
n — > n r and n — > mj., 



^ T(k/2,a 2 \n r + r| 2 ) ^ ira k 2 cos(27rra fc r) /7r 2 |ra fe | 2 \ 



,n=0 



(3.115) 



v T(k/2,a 2 \n r \ 2 ) v na k 2 ( 7r 2 |m fc | 2 \ _ a k 

n ^ r(fc/2)|n r |^ + J^ T{k/2) H a* j + %n " r( | + 1} 



(3.116) 



where ip^ is given by the expressions (3.71), (3.74) or (3.75). 

For a long-range interaction within the jellium model, the potential energy becomes 

U i6L = jiT l ^ iS /L ) + ^ (3.117) 



with the notation 



L = (L x L y ) l l 2 (3.118) 
n r = n x L x /L + n y L y /L (3.119) 
m k = m x L /L x + m y L /L y . (3.120) 



3.6. Ewald method for Yukawa potential 



71 



For a multicomponent gas (quasi-neutral in case of a long-range potential), the potential 
energy is 

^ SCn = Tfc£ m^a/Lo) + f-f £ . (3.121) 

Finally, the usual modification to calculate the momentum space sum linearly in N is 
given by 

a 2 1 ~ Na 2 ~ 

Uiel = JkE^/Lo) + — k £ x(m k )\S cq U™ k )\ 2 + -fft (3.122) 

L i<j ZL m k ^0 

U" cn = iZ^K/Lo) >c{m k )\S q {m k )\ 2 + f-f £ , (3.123) 

^0 i<j ^0 m k ^0 

with ijj, Scquai, S q defined by (3.109) - (3.112) in their corresponding two-dimensional 
variants. 

Some of non-orthogonal lattices can be sampled using the concept of pseudo-unit cell. 
For example a triangular lattice is constructed by translation of a single atom along two 
vectors with 60° angle between them. The same filling can be obtained by translation of 
a rectangular pseudoelementary cell with two atoms which can be readily calculated with 
the presented formulas. 

3.6 Ewald method for Yukawa potential 

As it was mentioned above, the Ewald summation technique can be applied to interaction 
potentials of a more generic kind, than polytropic ones, for example to the Yukawa class 
of interaction E p (ri — Vj = exp(— a\r\)/r. The interaction potential of this form widely 
appeared in nuclear physics as a primitive model potential inside nuclei and in simulations 
of plasmas, where it is used instead of the Coulomb potential to reflect the screening 
properties of plasma, and in the other applications. 

The derivation of Ewald sums in the style of de Leeuw et al. [dLPS80] , that we applied in 
3, cannot be used directly with the Yukawa term exp(—a\r\)/r, but its certain modification 
may actually be used. Let us briefly explain how it can be done, having in mind the general 
line of derivation, given in the chapters above. 

Consider the pairwise potential E p (r) = exp(— a\r\)/r, with a for a positive constant, 
defining a screening size. We can repeat all the procedure of obtaining the form of p(r) 
term (real-term component of the Ewald sum) unchanged for Coulomb potential with the 
coefficient exp(— a\r\). Therefore this term is obtained straightforwardly as 

p(ri) = exp(— a|n|)erfc(a|n|)/|n| (3.124) 

The other component £(m) cannot be found directly in this manner, since the Jacobi 
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transform, relying on a multidimensional variant of the Poisson formula [SW71] 

£/(r + n) =^F[/](m)e 2d ™ (3.125) 

n m 

with f(u) = exp(— u 2 — a\u\), requires for an analytical form of 3D Fourier transform 
F[f], which is unknown for a > 0. Nevertheless, it can be found computationally with an 
arbitrary precision, thus allowing to discover all the components of ip{ r ) an d £■ As the 
summation in the Fourier space is a summation of cosines with the factors x(m), depending 
only on the magnitude of the wave vector |m|, the latter factors can be precalculated and 
taken from a one-dimensional array. 

The other approach to the derivation of Ewald sums was presented in a number of 
works [SC00, Maz07a, Maz07b] and based on a traditional treatment of the Ewald tech- 
nique for Coulomb systems. The first step of the approach it to consider the charge 
distribution as 6(r) = (S(r) — d{r)) + d{r) = d\{r) + d2(r), where d{r) stands for a charge 
density, represented by a function, well localized around the point charge with suitable 
mathematical properties, normally it is a Gaussian. This background p(r) has to obey the 
condition, that the integral charge below the Gaussian is equal to the point charge, that 
is the background is neutralizing. The variance of the Gaussian is taken as an arbitrary 
parameter and is a subject of optimization in the following treatment. The charge density 
profiles are then treated separately, making use of the fact, that the Yukawa potential is a 
Green's function of a Helmholtz equation of a certain form. The final solution reduces to 
finding the Green's functions for two Helmholtz equations for charge densities d\(r) and 
g?2(V)> which is a relatively simple analytical problem. The periodic (Ewald) form Yukawa 
potential is then equal to the sum of these Green's functions. 

The final expressions of the Ewald sum for the Yukawa potential in our previous nota- 
tions {p(n), x(m), £, Ci, C2} can be written in the following form: 



erfc(a 


n 


+ a/(2a))e a W H 


- erfc(a 


n 
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p(n) = — — — — (3.126) 



= 7r eXP 2 l 2 a \ /A } (3.127) 
7r 2 m 2 + a 2 /4 

Cl = 4 7T eXp( ~^ ) (3.128) 
a 2 

2a ( a 2 \ „ / a 



C 2= exp -_ +.«fc^) (3.129) 



where a is dimensionless, when r is given in the units of L. Notice that the limit a — > 0, 
corresponding to a Coulomb system, yields incorrect results in the constants C\ and C2. 
This may be seen as a reflection of the fact that the initial sum over n for a Yukawa system 
with any finite parameter a is absolutely convergent, which is not the case for a system 
of Coulomb charges, thus here the continuous passage to the limit a — > is an incorrect 
operation. 



3.7. Summary of the analytic results 
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3.7 Summary of the analytic results 

In the previous sections, we have derived general expressions of the Ewald sums for poly- 
tropic l/|r| fe potentials in three- two- and one-dimensional systems. For integer values of 
k, the polytropic potential reduces to a power-law interaction, which comprises realizations 
of high physical relevance. Integer power-law potentials include 

• k — 1 - Coulomb l/\r\ interaction; 

• k = 2 - Calogero-Sutherland 1/|?"| 2 interaction; 

• k = 3 - isotropic 1/|t*| 3 component of dipole-dipole interaction (one-dimensional 
systems; two-dimensional system of dipole oriented perpendicularly to the plane); 

• k — 4, 5, 6 - interaction between different Rydberg atoms; 

• k = 6, 12 - van der Waals interaction. 

The expressions for the potential energy for both the jellium model and the general 
case of a charge-neutral simulation cell are the following 



n m^O 

£ = £ p{n r ) + + C x + C 2 

R(n, r) = p(n + r) 

K(n, r) = x(m) cos(27rmr) 

c m = ^ (fe _ 3)r[|] HA^d 

-47r + 4vrln(2a) if fc = 3 



C 2D = I (fc - 2) r[|] ifA; ^ 2 
27rln(a) ifjfe = 2 



C 2 
L 



a 



r(f + i) 



f (L x L y L z y/ 3 in 3D 
\ (LncLy) 1 / 2 in 2D 

n r = (n • L)/L , with X = (L x , Lj,, L z ) 

m fc = (m • witn If' = (l/L x , 1/L y , 1/L Z 



3.130) 

3.131) 
3.132) 
3.133) 

3.134) 
3.135) 

3.136) 

3.137) 
3.138) 

3.139) 

3.140) 
3.141) 
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Alternatively, by performing a momentum space sum the above set of equations become 

U scn = iSft^Cry/Lo) 4E K(rn k )\S q (rn k )\ 2 + |-f f (3-142) 
^0 i<j ^o mfe ^ Ziv o 

(7 2 1 ~ No 2 ~ 

U ' e = 7^£^/Lo) + £ x(m fc )|5 cqual (m fc )| 2 + -fe (3.143) 

^0 i<j ^0 m k ^0 £L 

${r) =Ei?(n r ,r) + C 1 (3.144) 

n 

£ =J2p( n r) + Ci + C 2 (3.145) 

^cquai = g- £exp(27rim fc r i /L) (3.146) 

j 

S g = ^g i exp(2vrim fc r i /L) . (3.147) 

3 

In accordance with considerations discussed in preceding sections, the simulation cell 
has to fulfil the charge neutrality condition (X)£Li Qi — 0) for long-range potentials. Also, 
notice that in the particular case of a cubic simulation cell, n r = n, m k = m. 

Explicit expressions of the coefficients p(n) and x(m) for the most relevant interactions 
are summarized for 3D and 2D systems in Table 3.1 and Table 3.2, respectively. 



Table 3.1: Coefficients p(n) and x(n) taken from Eqs (3.37) and (3.18) for 3D geometry. 
LR and SR stand for long range and short range, respectively. 
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The expressions of the Ewald summation components are consistent in appropriate 
limits with the earlier published results by M. Mazars [MazlO, Mazll] (inverse power-law 
potentials in 2D and 3D geometries), by C. Mora et al. [MPW07] (2D dipolar bosons), by 
N. Karasawa et al. [KI89] and W.-Z. Ou-Yang et al. [OYLS + 05] (Lennard- Jones potential). 



3.8. Practical application and optimizations in the Ewald technique 
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3.8 Practical application and optimizations in the Ewald 
technique 

3.8.1 Optimization scheme with N 2 dependence 

The basic idea of the Ewald method is to calculate slowly decaying sums in a rapid manner 
by means of the Fourier transform of the slowly converging part. Although conceptually it 
provides an exact result, the number of terms which has to be summed in order to reach 
the needed convergence is a priori unknown. Once we choose the interaction potential, 
this fixes the exact form of the sums to calculate, and the practical remaining question is 
the proper choice of the free parameter a and the numbers of terms to be calculated in 
the sums, originated from coordinate and momentum spaces. The total computer time T 
is obtained from the time needed to evaluate different sums 

T = (t r N r + t k N k )N 2 /2, (3.148) 

with the constants t r and t k depending on the complexity of the coefficients in the sums 
and the factor N 2 /2 approximating the number of pairs. Here N r and N k are numbers of 
terms which are summed for each pair r^, in particular the case N r = 1 corresponds to the 
so called minimum image convention. One can notice that t k is usually much less then t r , 
since in the Jacobi-transformed sum we only calculate cosine functions, which is generally 
far less time-consuming than the complicated functions appearing in R. It is clear that 
the parameter a affects only the resulting error in the energy. In fact, the value of a being 
very small or very large eliminates errors in one of the sums, but amplifies them in the 
other, so there is an "optimal" point for a, yielding a minimum error in the total energy. 



Table 3.2: The coefficients p(n) and x(m) taken from Eqs (3.37) and (3.77) for 2D geom- 
etry. LR and SR stand for long range and short range, respectively. 
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In the following, we discuss a way for error (SE) minimization assuming the calculation 
time T fixed. From our point of view, a useful approach for practical implementation is 
represented by the following scheme 

• We determine a time law T = (t r N r + t^N^N 2 /2 in a preliminary calculation and 
fix the values of t r and tk- 

• We take a set of configurations, corresponding to the equilibrated state using an 
initial Ewald summation. Then, we calculate the exact energies E ex (as a converged 
result of the Ewald summation) and the energies E(a, N r , N k ) biased by a choice of 
N r and iVfc. For each pair (N r , Nk), we find an optimal value of a = a opt (N r , Nk). 

• We choose the goal accuracy 5E acc (normally, well below the statistical error). We 
plot the error as a function of the computer time spent and choose the less time 
consumption case among the points that lie below SE &CC , therefore obtaining all the 
parameters required: a, N r and Nk- From now on, these parameters are used in 
actual simulations. 

3.8.2 Example of optimization 

Let us illustrate the scheme proposed in the preceding subsection taking as an exam- 
ple the problem of two-dimensional zero-temperature Bose gas of particles, interacting 
through the pairwise repulsive Cs/\r\ 3 potential. The model corresponds to the dipole- 
dipole interaction with all dipole moments aligned perpendicularly to the plane of mo- 
tion. The simulation is performed with N = 108 particles in a quadratic box (L x = 
L y = L = 1). The dimensionless Hamiltonian in the present example is taken to be 
H = — J2f=i Vf/2 + Y,i<j Vl r «il 3 w hh the unit of energy being h 2 /mL 2 . To describe the 
ground-state properties of the system we use the variational Monte Carlo (VMC) method 
and a Jastrow wave function with a two-body correlation factor which is solution of the 
two-body scattering problem [McM65]. 

The optimization is done by averaging over iV con f = 50 uncorrelated VMC configura- 
tions, sampled according to the chosen probability distribution. We define the error 5E(a) 
as a sum over iV con f configurations of the difference of the Ewald energy E(i coni) a, N r , N k ), 
calculated for a given set of parameters (a, N r , N k ) and the converged energy E cx (i con{ ) = 
lini/Vfc-Kx) hm^^oo E(i con f, a, N r , Nk). The dependence of the computer time T, needed for 
the evaluation of Ewald sums, on the parameter set is shown in Figs 3.1 and 3.2. In Fig. 3.1, 
we show the dependence of T on the number of terms N r in real space for different fixed 
numbers of terms Nk in the momentum space. The computation time is proportional to 
the number of terms and the resulting dependence is linear in N r . A fixed number of terms 
Nk requires a certain amount of calculations which results in a constant shift. Similarly, 
keeping N r fixed and varying Nk produces a linear dependence in Nk with a constant shift 
which depends on iVV, as shown in Fig. 3.2. 
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Figure 3.1: Dependence of the calculation time T on the number of terms N r in the 
coordinate space for fixed numbers of terms in the momentum space Nk = 5,9,25,45. 
Symbols are calculation times. Lines are linear fits for the data. 




Figure 3.2: Dependence of the calculation time T on the number of terms in the 
momentum space for fixed numbers of terms in the coordinate space N r = 1,5,9,21. 
Symbols are calculation times. Lines are linear fits for the data. 
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As one sees in Figs. 3.1 and 3.2, the time dependence is linear both on Nf. and N r , 
although the point corresponding to (0,0) in (N r , Nk) does not necessarily gives T = 0, 
since the reported time also contains some initializing calculations. The total error in the 
potential, as it is defined above, is given by 



According to our previous considerations, in the case of very small or very large values of 
a the error coming from one of two sums, that is in the real or momentum space, grows 
and dominates over the error coming from the other sum; for a certain "optimal" range 
of a these two errors are of the same order. Notice that for each particular configuration, 
and each pair (N r , N k ), it is possible to find a opt (i con f), such that E(a opt (i C onf),iconf) — 
E ex (iconf) = 0. Instead, our task is to obtain a "universal" parameter ao, minimizing the 
total error (3.149). The mean over the configuration set of the biased energies E(ato, i CO nf) 
is used as an estimation for the mean of the exact energies E ex , introducing an inevitable 
systematic error. As it appears in typical calculations, this error is at least one order of 
magnitude smaller than the statistical error (3.149) given by the minimization of 5E(a). 
In our benchmark calculations we also checked the dependence of the total energy on 
the value a for different pairs (N r ,N k ), which revealed characteristic plateaus for certain 
ranges of a (of the order of 1). It means that the Ewald summation indeed converges fast 
to a universal result. Nevertheless, the optimum value of the parameter a, minimizing the 
cumulative error, depends on the cut-off numbers in the sums in the range from a = 1.45 
(N r = 21, N k = 5) to a = 5.25 (N r = 1, N k = 45). 

A second step is the study of the dependence of the error and time on different pairs 
(N r , Nk). The calculation time can be split as the sum of times for summing up in 
coordinate and momentum spaces, T = (N r t r + N k t k )N 2 /2 as in Eq. (3.148), with N r , N k 
being the numbers of terms in each sum for a single pair. Every one of these sums converges 
when N r , Nk — > oo to a certain value, depending on a, while the sum of the limiting values 
is a constant. We can take into account the errors, corresponding to each of the sums 
separately. For a — > the error for the real space term is zero and the other one tends to 
infinity (and vice versa as a — > oo). The minimum total error should therefore correspond 
to the value of a, satisfying the relation d(5E r + 5E k )/ da = 0. 

Focusing on the 2D system of our example, we note that the long-range expansions 
of the terms in (3.37) and (3.77) are similar, in a sense that the leading terms in both 
expressions are Gaussians, 
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Figure 3.3: Resulting error of the energy per particle as a function of the computer time 
for different parameter sets. Symbols are the errors, related to truncation in the Ewald 
method. Lines are exponential fits for the data. 



3.8. Practical application and optimizations in the Ewald technique 
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The power-law terms in n, m and the constants C r , C& may be neglected since the leading 
behavior is driven by the Gaussian. The cut-off errors due to finite numbers of elements 
in the sums can be evaluated by ignoring the discrete structure of the images and approx- 
imating the sums by uniform integrals, 
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(3.152) 



with R ~ JN r /7r and K ~ ^N k /iT the approximate cut-off lengths in real and momentum 
spaces, respectively. The optimal value for a can be obtained by solving the equation 
dSE r / da = — d5E r / da. The first-order approximation of this equation is found by taking 
logarithms of both sides and omitting constants and terms, depending on a logarithmically, 
that is 

A k /a 2 - A r a 2 = (3.153) 



with A k = uN k and A r = N r /ir, which yields 

a = (A k /A r ) 1/4: = (n 2 N k /N r ) 
Then, at lowest order one finds (3.152), 



1/4 



N 2 „ N 2 i 

5E ~ — exp(-a 2 iV r /7r) = — ex P (-^fN k N r ) . 



(3.154) 



(3.155) 



Since the calculation time is linear with the numbers of elements N r and iV^, we may 
conclude that with N k fixed and comparatively large iV r , \n(5E) ~ ~ VT and vice 
versa, with N r fixed and large N k , hx(SE) ~ V^fe ~ VT. This power law may be easily 
checked in our calculations, as it is shown in Fig. (3.3). Note that for the obtained value of 
a the errors of the real- and momentum-space cutoffs are of the same order of magnitude, 
that is SE r w 5E k , which may serve as a rough criterion to optimize the parameter a. 



3.8.3 Optimization scheme with TV 3 / 2 dependence 

A more advanced procedure for optimization of the parameters, proposed by Perram et 
al. [PPdL88], yields an asymptotic scaling iV 3//2 , with N the number of particles. It is 
based on the form of Ewald summation with the rearranged momentum space sum, linear 
in N (3.142). Note that here the linear N dependence is obtained after a summation over 
particle coordinates, while in Sections (3.8.1) and (3.8.2) the momentum space sum was 
evaluated over the pairs of the particles and had N 2 dependence. Suppose the values of the 
calculation time t r , t k to perform unit computations in both sums are known and the target 
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error level exp(— p) is fixed. Then, the total execution time in the real and momentum 
spaces is 

T = T r + T k = N\R\ + NirK% (3.156) 

with p = a 2 R 2 = 7r 2 K 2 /a 2 . Expressing K as K = p/(irR) we can see that the minimum 
of the total time T corresponds to 

"opt = 



^) N 1 ^ . (3.159) 



The computation time is equally divided between the real and momentum space parts (this 
was also stated in our simple optimization scheme), with a scaling of the whole summation 
given by 

T = 2N 2 7iR% = 2p v / WV 3/2 (3.160) 

Notice that the values of the free parameters change very slowly when the simulation cell is 
enlarged, and in particular a is not affected by the choice of the precision. Similar formulas 
for the optimized parameters in three-dimensional systems, with a discussion of different 
techniques to improve performance of the Ewald summation, are given by Fincham [Fin94]. 
A more precise and detailed analytic study of the cut-off errors with verifications of the an- 
alytic results in actual calculations can be found in the work of Kolafa and Perram [KP92]. 
An optimized method for treating the truncation error in Ewald sums with generic po- 
tentials was proposed by Natoli and Ceperley [NC95]. While the needed CPU time scales 
as 0(N In iV) 3 / 2 , it was shown that in the example of the Coulomb potential the method 
resulted in greatly improved accuracy compared to that of standard Ewald technique for 
a comparable computational effort. This method is based on an expansion of the real 
space function in an arbitrary radial basis with a parametric set of numbers in place of 
the fc-dependent prefactors of exp(27rimr). The subsequent minimization of \ 2 with re- 
spect to the whole set of parameters yields a final optimal solution, that is the real space 
expansion coefficients and the fc-space factors. This technique was also applied to derive 
the optimized summation formulas for the two-dimensional Coulomb system [HB05]. 

In general, the unit computation time in momentum space is 2-4 times faster than the 
one in real space. Taking the following reasonable assumptions p = 4tt, tk/t T = 3, we 
find R opt « 2.6/iV 1 / 4 (the box size is taken to be 1). We want R to be below 0.5, since 
in this case the summation in the real space reduces to the accumulation of the single 
component n = ("minimum image convention"). This condition R opt = 0.5, with our 
previous assumptions, corresponds to 

iV opt = 770, fT opt = 8.0, o opt = 7.1 (3.161) 



3.9. Conclusions 
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In smaller systems, the other components of the real sum, starting from \n\ = 1, should 
be considered. 

It is worth pointing out that if the interaction is very strong at short distances (as 
for the Lennard- Jones potential), then in principle the real-space cut-off R can be chosen 
below the "hard core radius" with a large enough value of a. This leads to the possibility 
of dropping completely the real-space part of the total sum and treat the /c-space only 
This can be advantageous in different aspects, especially with the current progress in the 
development of efficient FFT-based methods. Nonetheless, we are not aware of any present 
application of a similar technique. 

3.9 Conclusions 

In the present Chapter, we have applied the Ewald summation method to l/|T"| fc polytropic 
potentials in three-, two- and one-dimensional geometries in a simulation box with periodic 
boundary conditions. We have found the explicit functional forms for all the components 
of the sums in both real and momentum spaces, with special attention being paid to the 
cases of long-range interactions, that is conditionally convergent or divergent potentials 
(corresponding to k < D, with D standing for the dimensionality), "marginal" interactions 
(k = D), and short-range interactions (with k > D). For the latter case of short-range 
interaction potentials, where in principle a straightforward summation of the initial sum 
(3.1) is possible, the Ewald method is shown to be useful, as it yields the faster (Gaussian) 
convergence rate. A condition of charge neutrality of the simulation cell is stated to be 
necessary for conditionally convergent and divergent potentials; a homogeneous positive 
charge background ("jellium" model) is introduced as the most relevant and frequently 
used kind of neutralization. The conditionality of the convergence for a charge-neutral 
system, governed by the Coulomb interaction, is discussed with a justification of the use of 
a specific periodicity-preserving convergence factor. The derivation technique, presented 
in our work, is consistent with the arguments of de Leeuw et al. [dLPS80]. 

The results are first presented for the case of a 3D system in a cubic simulation box 
in order to explain the general mathematical procedure, which for the specific case of the 
Coulomb potential recovers well-known results [AT89]. Later on, the same mathematical 
technique is applied to 2D and ID geometries. For the one- dimensional case the initial sum 
for the potential energy is explicitly evaluated (3.96), nonetheless the Ewald summation 
is developed for this case too and may be used as a mathematical equality. The special 
representations of the reciprocal space sums, linear in the number of particles N and hence 
more efficient in actual modeling, are presented for 3D and 2D systems. The explicit 
expressions for the terms of the Ewald sums are given in a tabular form for physically rel- 
evant potentials with small integer power indexes k, as dipole-dipole interaction potential, 
Lennard- Jones potential and others in both three- and two-dimensional geometries (see 
Tables 1 and 2). 
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When the simulation box cannot be chosen cubic, for example in a modeling of a 
three-dimensional hep crystal structure, the Ewald method can also be applicable after 
a certain modifications. Formally, it consists in the choice of an appropriate rectangular 
simulation box and a substitution of the vector n by n r = (n x L x + n y L y + n z L z )/L Q and 
rrik = (m x /L x + m y /L y + m z /L z )L in the real and momentum space sums, respectively 
[see (3.102) and (3.103)]. 

The optimization of the involved parameters, that is cut-off numbers in both sums and 
the integration parameter a, is a necessary operation in order to improve the convergence 
rates and avoid excessive calculations. The main idea of the optimization, proposed in 
the present work, is to perform a benchmark calculation, minimizing the variance of the 
result. A particular example of the application of the technique is presented for a calcula- 
tion of the potential energy of a two-dimensional gas of dipoles, aligned perpendicular to 
the plane of motion. This practical optimization technique is thought to be efficient for 
stationary and nearly uniform systems that appear, for instance, in Monte Carlo simula- 
tions. In spite of being very simple, it allows to find rather quickly adequate parameter 
ranges. The analytical estimations of the parameters are given as well and are proven to 
be consistent with the results, obtained in our method. A more sophisticated method to 
optimize the calculation parameters, taking advantage of the O(N) representation of the 
Fourier transform sum, is also presented with explicit estimations of the parameters for a 
typical system simulated by Quantum Monte Carlo methods. 



Chapter 4 

Phase diagram of a Yukawa system 



4.1 Introduction 

Recent advances in trapping and controlling ultracold dilute gases have permitted to realize 
highly tunable and extremely pure Fermi systems [DGPS99, GPS08]. This has provided 
new insight in the study of fundamental problems in condensed matter physics. For ex- 
ample, the original BCS theory [BCS57] was developed to explain superconductivity in 
metals, where the control over interactions and densities is very limited. However, in re- 
cent experiments with ultracold Fermi gases in the BCS-BEC crossover the strength of 
the interactions is controlled by external magnetic fields in the vicinity of a Feshbach res- 
onance, while the geometry is tuned by means of magnetic or optical confinement. This 
has allowed, for instance, to measure the equation of state in the BCS-BEC crossover 
in high precision experiments [NNCS10, KSCZ12]. Numerically, the best calculation of 
the zero-temperature equation of state is obtained in quantum Monte Carlo simulations 
[CCPS03, CPCS04, ABCG04, FGG11, CGSZ11] 

After the big success achieved with single species there is nowadays a growing inter- 
est in fermionic mixtures. Quite recently, fermionic mixtures consisting of atoms with 
different masses have been realized experimentally [TVA + 08, IKH + 11] and studied theo- 
retically [GGSC09, BRS09, BD10]. Novel physical phenomena like Efimov states [Pet 03, 
HHP10, LP11, YZZ11, WLvSE12], trimer and cluster formation might be observed [KMW+06, 
KFM+09, ZDD+09, BWR+09, WHH+09, WLO+09, KFB+09, NHM+10] in these systems. 
The case of large mass imbalance is especially interesting, and mixtures of 6 Li and 40 K 
are being investigated experimentally [TVA+08, WSK+08, VTC+09, STN+09, TGL+10, 
CBV + 10, TKZ + 11, RCS + 11]. Even larger mass ratios are reached in mixtures of 6 Li and 
173 Yb [IKH + 11, HTY + 11]. In this chapter we present results for the phase diagram of 
Fermi mixtures as a function of the mass ratio using quantum Monte Carlo methods and 
determine how crystallization of this system can be realized. 

From the theoretical point of view, it was proposed in Ref. [PAP + 07] that an effective 
Yukawa interaction, induced between heavy-light pairs of fermions, might lead to crys- 
tallization in quasi-two-dimensional systems. In this work we extend that discussion and 
analyze the possibility of realizing a gas-crystal phase transition at zero temperature in 
three-dimensional systems. We obtain the phase diagram and discuss how large mass ratios 
have to be for reaching crystallization. 

The interest in the phase diagram of quantum Yukawa particles is rather old as the 
Yukawa potential has long been used, for instance, as a model for neutron matter [BP75, 
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Kaw81]. The Yukawa potential also describes interactions in dusty plasmas where charged 
dust particles are surrounded by plasma which introduces screening [HFD96, HFD97, 
SAll]. The Yukawa potential is often used as well as a model for suspensions of charged 
colloidal particles [KRG86, RKG88, RT87, MF91, HHBN11, GNP12]. The classical finite 
temperature phase diagram has been extensively studied [HFD96, HFD97, SAll, KRG86, 
RKG88, RT87, MF91, HHBN11, GNP12] while much less is known about the full quantum 
phase diagram. 

In the 70's, Ceperley and collaborators [CCK76, CCK78] used the diffusion Monte 
Carlo algorithm to estimate the zero-temperature phase diagram of the Yukawa Bose fluid. 
In their work the phase diagram was built assuming that the Lindemann ratio remains 
constant along the solid-gas coexistence curve, with the explicit value being evaluated only 
at a single point. In the present Chapter we carry out a full study of the transition curve 
and present the phase diagram in terms of experimentally relevant densities and mass 
ratios of heavy to light fermions. The Lindemann criterion prediction has turned out to 
be quite precise apart from the region of high densities. 



Mixtures of fermions with different masses have been realized recently in a new genera- 
tion of experiments [TVA+08, WSK+08, VTC+09, STN+09, TGL+10, CBV+10, TKZ+11, 
RCS + 11, IKH + 11, HTY + 11]. The interactions can be tuned to allow the formation of two- 
component molecules. The s-wave interactions within a single component are prohibited 
due to Pauli principle. Yet, an effective interaction between same-spin fermions can be 
induced by the presence of the other component. The limit of large mass ratio has been 
analytically addressed in Ref. [PAP + 07]. The effective interaction between heavy particles, 
which was obtained in the limit of large distances within first Born approximation, has the 
form of a screened Coulomb (Yukawa) potential. This leads to a description of the system 
in terms of a composite (molecular) bosonic gas interacting with an effective potential. 

We study a system of heavy fermions of mass M interacting among themselves and 
moving on a background of light fermions of mass m. The net effect induced by the 
movement of the light fermions can be characterized by a Yukawa potential, leading to 
the following effective Hamiltonian [PAP + 07] describing the interaction between composite 
bosons formed by pairs of heavy and light atoms 



where a is the atom-atom s-wave scattering length between two atomic species and r^ 
are positions of heavy atoms while the positions of light atoms have been integrated out. 
The ground-state properties of the system are then governed by two dimensionless pa- 
rameters, namely the gas parameter na 3 and the mass ratio M/m. Equivalently, Hamil- 
tonian (4.1) describes a bosonic system interacting via the screened Coulomb potential 



4.2 Model Hamiltonian 
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Vi n t(r) = gexp(— Xr)/r by mapping the charge to q = 2h 2 /ma and the screening length to 
A" 1 = a/2. 

We calculate the ground-state properties corresponding to the Hamiltonian (4.1) by 
means of the diffusion Monte Carlo (DMC) algorithm [BCN94]. This method solves 
stochastically the Schrodinger equation in imaginary time providing the exact energy 
within controllable statistical errors. The coexistence curves can then be traced by di- 
rect comparison of the energies of the solid and gas phases. The efficiency of the DMC 
method is greatly enhanced when importance sampling is used. This is done by multiply- 
ing the (unknown) ground-state wave function ip(ri, . . . , rjv) by a guiding wave function 
iPt(ti, • • • , rjy) and solving the equivalent Schrodinger equation for the product. As a result, 
the points in phase space where the guiding function is large get sampled more frequently 
and this improves convergence to the ground state (see Section 2 of this Thesis). 

The properties of the gas phase are studied by constructing the guiding function in 
a Jastrow two-body product form i/)t(ti, ■ ■ ■ , rjv) = Tii<j hd^i — ry|). We determine 
the optimal two-body Jastrow term /2( r ) by solving the corresponding Euler-Lagrange 
hypernetted-chain equations [Kro98] (HNC/EL), discussed previously in Section 2.4.3.4, 
discarding the contribution of the elementary diagrams. The resulting wave function cap- 
tures basic ingredients coming both from the two- and many-body physics of the problem. 
On the other hand, the energy of the solid phase is obtained by using a Nosanow- Jastrow 
guiding wave function tpri^i, ■ ■ ■ , ?n) = UiLi fi( r i — r- att -) Yl i<: j /2(l r i — r il) with Gaussian 
one-body terms /i(rj — r- att ') = exp(— a(r — r' att ) 2 ) describing the localization of particles 
close to the lattice sites *. The parameter a controls the localization strength and is 
optimized by minimizing the variational energy. 

In order to find the energy in the thermodynamic limit, we carry out simulations of a 
system of N particles in a box with periodic boundary conditions, and take the limit iV — > 
oo while keeping the density fixed. In the simulation of the crystal the number of particles 
should be commensurate with the box which restricts the allowed number of particles. 
For fee packing the simulation box supports N = 4z 3 = 4, 32, 108, 256, .... In order to 
add more values we also use periodic boundary conditions on a truncated octahedron (see 
Appendix B), which allows simulations with N = 2i 3 = 2, 16, 54, 128, 250, 432, . . . particles 
with a larger effective volume of the simulation box and reduced anisotropy effects. Finally, 
the convergence is further improved by the Ewald summation technique [Ewa21, OAB12] 
in the cubic box, which we use in the calculations at large densities. 

In Fig. 4.1 we show two characteristic examples of the finite-size dependence of the 
energy at two different densities. For large enough system sizes, the energy is well fitted 
by a linear dependence in 1/N. For small number of particles the behavior is no longer 
linear, especially at large densities due to strong interparticle correlations. We find that 
system sizes of N > 100 have to be used in order to ensure the linear regime at considered 
densities. The thermodynamic energy is then obtained as a result of a linear extrapolation 
1/N ->■ 0. 



88 



Chapter 4. Phase diagram of a Yukawa system 



0.999 - 




0.998 - 



0.997 



n a = 1.6 



♦ n a = 0.192 



0.000 



0.005 



UN 



0.010 



Figure 4.1: An example of finite-size dependence of the energy in the gas with periodic 
boundary conditions in truncated octahedron for M/m = 187 at two different densities 
na 3 = 1.6 (upper set of data points) and na 3 = 0.192 (lower set of data points). Symbols, 
DMC energy; lines, linear fit to energy for large system sizes. Energies are scaled with the 
thermodynamic value E extr , obtained in 1/N — > extrapolation. 
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4.3 Phase diagram 

An intrinsic property of Coulomb particles is to self-assemble into a Wigner crystal at low 
densities and to remain in a gaseous phase in the opposite limit, due to the long-range 
character of the interaction [Wig34]. The Yukawa potential is similar to the Coulomb 
one at densities large enough for the interparticle distance to be much smaller than the 
screening length, which is fixed by the s-wave scattering length a between two different 
species of atoms. One then concludes that the Yukawa system stays in a gaseous phase 
at large densities. In the opposite regime of small densities, na 3 <C 1, the interaction 
potential decays exponentially fast showing a short-range behavior that leads the system 
to a gaseous phase. For example, the fee crystal of hard-sphere bosons of diameter a s melts 
at density na 3 s pa 0.24 [HLS71, KLV74, DNRA90]. The intermediate regime na 3 « 1 is 
the most interesting one, as crystallization may or may not take place depending on the 
strength of the interaction, which in the current case of the Hamiltonian in Eq. (4.1) is 
governed by the mass ratio M/m. A relevant question then is what is the minimal mass 
ratio at which crystallization can be observed. 

In order to obtain an accurate description of the phase diagram, we study the finite 
size dependence and extrapolate the energy to the thermodynamic limit. As we mentioned 
in the sections, devoted to the methodology, the diffusion Monte Carlo method is also 
generally biased by the average size of the population of walkers N w (sets of N p particle 
configurations), that can be controlled in the simulation. In practice the results of different 
calculations converge to a stable value, when N w is large enough, although the convergence 
can be achieved with much less N w with a better trial wave. Usually the population of 
walkers, when the convergence is reached, is between 250 and 500 and is not affected much 
by the size of a system, therefore a preliminary simulation can be carried out fast with 
a small number of particles. Figure (4.2) demonstrates the analysis of N w convergence 
for the Yukawa system of 64 particles in a liquid phase (no size correction added). The 
convergence is seen to be reached for N w ~ 200. 

Throughout all the DMC simulations that we run we use the second-order in time step 
approximation of the Green's function of the Hamiltonian, which means that the time step 
bias in the results is also 0(t 2 ). This dependence can be observed by performing a set of 
trial calculations in order to find a value of At, when the time-related error complies with 
our accuracy goal. In Figure (4.3) we plotted the energy per particle versus time step and 
extrapolated the data to by a parabola. Here we can suggest that the acceptable time 
step is approximately 1000 (the required accuracy level is around the statistical noise for 
each single calculation). 

The resulting energies of the gas and solid phases are then analyzed using the double- 
tangent Maxwell construction which provides the melting and freezing densities. The 
zero-temperature phase diagram parameterized in terms of the dimensionless density na 3 
and the mass ratio, is shown in Fig. 4.4. We find that for mass ratios smaller than the 
critical value M/m pa 180 the gas phase is energetically preferable at any density. On the 
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Figure 4.2: Dependence of the energy per particle E/N p on the inverse number of walkers 
1/N W for the Yukawa system in a liquid phase, 64 particles. 



other hand, for larger mass ratios there is always a gas-solid transition at low densities and 
a solid-gas transition at large ones. Energetically, both the fee and bec lattices are possible 
in the solid phase. It is very difficult to discern numerically which packing is preferred as 
the energies in different crystalline phases are extremely close. Still, in the large potential 
energy limit, corresponding to a mass ratio M/m 3> 1, it is enough to compare the potential 
energy of the classical crystals with different packings. A simple, geometrical construction 
assuming that particles are tightly tied to their equilibrium positions leads to a transition 
density na 3 m 1.58. This prediction is depicted as a blue dashed line in Fig. 4.4. In the 
low-density limit we numerically find the value of the s-wave scattering length a s of the 
Yukawa potential (4.1) and fit it as a s /a = 0.436 ln(M/m) with accuracy below 1% in the 
region of interest. Note that a is the s-wave scattering length of fermionic particles which 
lead to the effective bosonic Hamiltonian (4.1) while a s is the s-wave scattering length 
between bosonic Yukawa particles. For the sake of comparison we also plot in Fig. 4.4 the 
gas-solid transition line of hard spheres of size a s given by M/m = exp(1.424/(n 1//3 a)). 

The figure also shows the results of Ceperley et al. [CCK76, CCK78] which were ob- 
tained by doing DMC calculations for three characteristic points in the phase diagram close 
to the solid-gas transition line. Overall, the agreement between that prediction and our 
results is good, the main differences affecting the region of large density where Coulomb ef- 
fects are strong. To our best knowledge this is the first time that the high-density quantum 



4.4. Large mass ratios 



91 




Ax 

Figure 4.3: Energy per particle E/N p (green points with errorbars) versus the time step 
At for the Yukawa system in a liquid phase, 64 particles; red solid line for a quadratic fit. 
Energy units are E\ = h 2 /ma 2 , time units are h/Ei ■ 10 -6 . 

solid-gas phase transition is observed in a simulation of Yukawa systems. 

In the case of the fermionic molecules, the resulting critical mass ratio is much larger 
than Mj m pa 13.6 for which the system is unstable due to formation of Efimov states [Pet 03, 
HHP10, LP11, YZZ11, WLvSE12]. The obtained phase diagram describes properties of 
metastable fermionic molecules while the true ground state corresponds to a many-body 
bound state. The stronger the effective interaction is (that is, the larger the mass ratio), 
the more distant are heavy fermions and the smaller the overlap with localized Efimov 
states is. 

4.4 Large mass ratios 

According to our results, the minimal mass ratio for which the crystalline phase can exist 
is M/m ~ 180 and it is achieved at the somewhat large value of the gas parameter na 3 ~ 
0.3. At these densities the fermionic nature of the molecules becomes important as the 
Hamiltonian (2.65) is derived under the assumption that na 3 < 1/8 [PAP + 07]. Our bosonic 
model is expected to be reliable at smaller densities where the critical mass ratio is further 
increased. 



92 



Chapter 4. Phase diagram of a Yukawa system 



400 



5 



300 - 



100 - 



i i i 1 1 1 1 1 1 



i i i i 1 1 1 1 



i i i i 1 1 1 1 1 i 



200 - i 



gas 





0.01 



— •— DMC 

Ceperley (DMC + Harmonic) 

FCC / BCC classical crystals 

hard sphere 

• ««.»»««!" ......... I, , | . , , . ...... .1 



0.1 



1 



10 



n a 



Figure 4.4: Zero-temperature phase diagram of the Yukawa potential corresponding to 
the Hamiltonian in Eq. (4.1) in terms of the gas parameter na 3 and the mass ratio M/m. 
Red symbols: transition point as obtained from the double-tangent Maxwell construction 
applied to the Monte Carlo data energies extrapolated to the thermodynamic limit; dashed 
line: critical density na 3 = 1.58. . . at which the energy of perfect fee and bec packings 
are equal; dash-dotted line: prediction of Ceperley et al. [CCK76, CCK78] obtained by 



imposing a constant Lindemann ratio; short-dashed line: 
DNRA90]. 



0.24 [HLS71, KLV74, 
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The mixtures of different fermionic atoms have already been successfully realized in 
experiments [IKH + 11, HTY + 11] but at significantly smaller mass ratios. Probably, the 
largest directly achievable mass ratio currently is that of Yb and Li atoms, M/m = 29, 
which is still much smaller than the critical mass ratio needed to observe the formation of 
an ultracold crystal. 

An alternative way to realize a fermionic mixture with a large and variable mass ratio is 
to confine one of the components to an optical lattice. At low filling fraction the distances 
between atoms are large compared with the lattice spacing, and the separation of length 
scales allows the description of the movement of a particle in the lattice as that of a 
quasiparticle with an effective mass moving in a medium where the lattice is absent. In 
a deep lattice interactions between particles are much weaker than the confining energy 
and so, to a first approximation, one can consider that as the problem of a single particle 
diffusing in the lattice. 

An optical lattice created by counter-propagating laser beams imposes an external 
potential V^tt-O^) y, z) = V (sin 2 kx + sin 2 ky + sin 2 kz) on every particle. The diffusion of 
a particle over a large distance is then governed by the tunneling rate between neighboring 
sites. The diffusion is largely suppressed (and the effective mass greatly increased) when 
the amplitude of the optical lattice is large, i.e. when Vq ^> E r with E r = h 2 k 2 /2m the 
recoil energy. The excitation spectrum in the lowest band can be described by Bloch waves 
of quasi-momentum q and energy e (q) = §ftw — 2 J (cos q x d + cos q y d + cos q z d) + . . . with 
d = n/k the lattice constant [BDZ08]. At small momenta the spectrum is quadratic in q 
and can be interpreted as the spectrum Eo(q) = Eq + h 2 q 2 /2m* of a free quasiparticle with 
an effective mass m*. Within the lowest band approximation the effective mass is inversely 
proportional to the hopping parameter J, 



The tunneling is greatly suppressed in the deep optical lattice limit Vq ^> E r . To better 
understand the contribution of the tunneling term in the present case, a semiclassical treat- 
ment within the Wentzel-Kramers-Brillouin(WKB) approximation can be used to calculate 
the tunneling probability p. One finds that it is proportional to 



where X\ and X2 are the classical turning points. In the deep optical lattice limit one can 
assume V(r) —Ek, V(r), with x\ and X2 corresponding to the positions of two neighboring 
minima. The result ing integral can be easily evaluated and predicts an exponential form 
J oc exp(— JVqJEt). A more precise expression can be obtained from the width of the 
lowest band in the ID Mathieu-equation [BDZ08], yielding 




(4.3) 
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This expression, together with Eq. (4.2), provides an analytic approximation for the effec- 
tive mass m*. 

In order to determine the dependence of m* on the lattice parameters in a non- 
perturbative way we evaluate the diffusion constant D of a real particle moving on the 
lattice and compare it to the diffusion constant Do = h 2 /2m* of a free quasiparticle 
of effective mass m*. The diffusion constant is obtained by means of DMC propaga- 
tion in imaginary time by measuring the mean-square displacement ((r(r) — r(0)) 2 ) = 
((x(r) - x(0)) 2 ) + ((y(r) - y(0)) 2 ) + ((z(r) - z(0)) 2 ) where r = (x,y,z) denote particle 
coordinates. The diffusion constant is then extracted as D = lim h((r(r) — r(0)) 2 )/(2r(i), 
where d = 3 is the system dimensionality. The resulting dependence of m* on the lattice 
amplitude is shown in Fig. 4.5. The figure shows the Monte Carlo prediction (solid line) 
compared with the approximation of Eq. (4.2) with J taken from Ref. [BDZ08] (circles) and 
from Eq. (4.4) (dashed line). As it can be seen, there is an almost constant shift between 
m* obtained in the Monte Carlo simulation and Ref. [BDZ08] compared to Eqs. (4.2-4.4). 
We have found that the description in the relevant region of interest is very much improved 
by subtracting a constant shift E^ = —?>/AE r from Vq in the argument of Eq. (4.4). This 
last prediction is shown by a thin line in Fig. 4.5 and provides a good approximation for 
V > 10£ r . 

One can understand these results in the following way: in the absence of the optical 
lattice the effective mass and the bare mass coincide, so m* = m. As the amplitude Vq of the 
lattice is increased, the particle movement is slowed down and the effective ma ss incre ases. 
In the deep optical lattice limit the effective mass grows as m*/m oc exp(y / Vo/-E' r ) and 
so the ratio can be made arbitrarily large by increasing the amplitude Vq (for instance 
m*/m ~ 1000 at Vq/E t = 40; see the inset in Fig. 4.5). This mechanism allows for 
increasing the mass of one of the two components while keeping the other one unaltered, 
so that the ratio M/m of the fermionic mixture can be made as large as desired when the 
mass of the heavy component is identified with the effective mass m* . Consequently, and 
according to the phase diagram shown in Fig. 4.4, there is a wide range of densities where 
one could find the system in the crystalline superlattice phase. Heights of optical lattices as 
large as (35 — GO)^ are readily achieved in current experiments [WTL + 06, TWL + 06] and 
correspond to sufficiently large effective mass ratios for the crystallization to be realized. 

Both small density and large density transition lines are accessible for Yukawa interac- 
tion caused by screening in dusty plasma, colloids and neutron matter. On the contrary, in 
two-component Fermi gas only the left part of the phase diagram can be realized since the 
effective Yukawa interaction is valid only at low densities. In fact, the validity criterion for 
the interaction potential in Eq. (2.65) was studied in Ref. [PAP + 07] and was found to be 
well satisfied for distances larger then r « 2a which leads to the condition pa 3 < 1/8 when 
r is identified with the mean interparticle distance. In this way, for example, for pa 3 = 0.1 
and mass ratio M/m = 300 the system is expected to be in a crystalline form. Much larger 
effective mass ratios can be achieved for realistic [WTL + 06, TWL + 06] lattice heights of 
(35 — 60)E r . We thus conclude that by using an optical lattice, a fermionic mixture of 
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Figure 4.5: Effective mass as a function of the lattice amplitude Vq in units of the recoil 
energy E r . Solid line: results obtained from the diffusion constant evaluated by propagation 
in imaginary-time; circles: lowest band approximation of Eq. (4.2) with values of J taken 
from [BDZ08]; dashed line: same results with J from the expansion in Eq. (4.4); dash- 
dotted line, same expansion with Vo shifted by — 3/AE r . Inset: same results on a semi- 
logarithmic scale. 



very different mass components can be used to test the phase diagram of the equivalent 
Yukawa model. 

4.5 Conclusions 

To summarize, in this Chapter we have obtained the zero-temperature phase diagram 
of bosons interacting through Yukawa forces. We have used a diffusion Monte Carlo 
simulation starting from a very good approximation to the optimal variational ground- 
state wave function obtained by solving the corresponding Euler-Lagrange hypernetted 
chain equations. The resulting phase diagram is very similar to the one originally obtained 
by Ceperley and collaborators [CCK76, CCK78], although significant differences arise at 
large densities. The phase diagram shows that any fermionic mixture of pure elements will 
always be seen in gaseous form, as the mass ratios required for crystallization of weakly 
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bound fermionic molecules are far beyond the ones that can be achieved in nature. Finally, 
we investigate an alternative mechanism based on the confinement of one of the species to 
a deep optical lattice which exponentially increases its effective mass as a function of the 
confining amplitude. The resulting mass ratio of the mixture created in this way can then 
be tuned at will and could be used to check experimentally the predicted phase diagram 
both in the gas and crystal (superlattice) phases. 



Chapter 5 

Phase diagram of Rydberg atoms 



5.1 Introduction 

Rydberg atoms have one electron excited to a high energy level. Such atoms exhibit 
strong and highly tunable interactions which may have an extraordinarily long range. Op- 
tically excited from suspended clouds of cold atoms, Rydberg atoms interact both between 
themselves and with the surrounding unexcited atoms, resulting in a rich behavior of the 
Rydberg systems. 

Due to the strong interactions, a Rydberg atom shifts the levels of nearby atoms suf- 
ficiently to prevent their subsequent excitation. A large number of studies deal with a 
local blockade regime. In such a regime a Rydberg atom blocks excitations in its vicinity, 
and the atomic clouds may be injected with well over 10 3 Rydberg excitations before the 
existing excitations block any further ones [HRB + 07, TFS + 04, SRLA + 04]. Unfortunately, 
the arrangement of the excited atoms in such experiments is not directly accessible and 
has been a subject of intense investigation. Understanding the ordering of Rydberg atoms 
may be important for interpretation of the experimental results, for example for the an- 
tiblockade effect predicted in [APPR07]. It was also suggested that a spatially ordered 
state may allow for a better control over quantum states in such experiments [PDL10]. 
Finally, there is an exciting possibility of observing phase transitions in these versatile 
systems, especially to states with long-range ordering [WLPB08, LWK+09]. 

Quantum many-body treatments attempting modelling of realistic Rydberg systems 
have been developed in the past [TFS+04, RH05, APPR07, SC10, PDL10, YRP+09], and 
were successful in reproducing a number of important experimental features [TFS + 04, 
YRP+09, LWK+09, AGHW10, SGH+10]. Due to complexity, it is often difficult to consider 
long-range order with such calculations. Nonetheless, strong short-range spatial correla- 
tions between Rydberg atoms were obtained in the calculations of Refs. [RH05, AGHW10, 
SC10], as the atoms avoid each other due to the blockade. Successful observation of 
the antiblockade effect was also a demonstration of a creation of the strong short-range 
correlations [AGHW10]. Possibility of long-range ordering (crystallization) of Rydberg 
atoms was recently predicted for systems coupled to specially selected chirped laser pulses 
[PDL10, vBSvL+11]. Ordering was also considered, and crystalline phase found, in theoret- 
ical calculations of both one and two-dimensional optical lattices [SLMD10, WB10, JAL11, 
SPG 11]. Remarkable non-commensurate crystalline phases in optical lattices emerged 
in Ref. [WB10]. 



98 



Chapter 5. Phase diagram of Rydberg atoms 



Given the complex nature of the interactions in the Rydberg systems, it is important to 
know how much of the behavior of large assemblies of Rydberg atoms stems directly from 
the pair potential of the interaction between the atoms. For this reason we aim to study 
ordering in the simplest model of the Rydberg systems. Because of the large number of 
Rydberg-excited atoms in the experiments, we consider the thermodynamic limit. While 
the results are established in the thermodynamic equilibrium, many present experiments 
with Rydberg atoms are too short to reach equilibrium. Thus comparison in such cases 
must be made cautiously. 



5.2 Model and methods 

The dominant interactions in the Rydberg systems are usually the Forster-resonant dipole- 
dipole interactions between the excited atoms. It was shown by Walker and Saffman 
[WS05, WS08] that, given a pair of Rydberg atoms in the same state, the interaction will 
not have zeroes as a result of the hyperfine structure or alignment of the atoms only if the 
resonant coupling is from the s to p states. Furthermore, interactions in the s + s — » p + p 
channels depend only weakly on the hyperfine structure of the p states, resulting in a 
nearly isotropic interaction, to within 10~ 2 . This perhaps in part motivates the use of 
the ns Rydberg states in current experiments [HRB+07, HRB+08, LWK+09, SGH+10]. 
Neglecting the hyperfine structure, the interaction for this resonance is isotropic and its 
matrix element is given in terms of the Forster defect 8 as [WS05] 



V(r) = S --si g n(6) ] J(^) + (5.1) 

which changes from V = C3/Y 3 to van der Waals' V = C^/r 6 (with Cq = —C$/5) for 
distances much larger than the crossover R^e ~ (— Cq/5) 1 ^ 6 . In the case of a strong local 
blockade, the blockade radius is often larger than the crossover distance. In such a case, 
the excited atoms are more likely to be found at distances where the interaction is already 
of the van der Waals type. 

The above arguments motivate the repulsive van der Waals model for the Rydberg 
atoms in the local blockade regime. We disregard any energy transfer or interactions with 
the underlying gas of the ground-state atoms, and particles are treated as spinless bosons 
in three-dimensional space with the many-body Hamiltonian 



i<3 I 7 "* r il 

Defining the reduced units of length and energy as 



ro ^_j .E^—^ (5.3) 
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allows to describe the properties of this model universally in terms of just two parameters, 
the dimensionless density prjj and temperature k^T / E . The units are selected to satisfy 
Eq = h 2 jmr\ = Cq/tq. The mass m in Eq. (5.3) is the mass of the atom. 

It is important to establish the applicability of the bulk phase diagram to finite systems. 
For a cloud of size R and number density p, the tail potential energy per particle can be 
estimated as pC§jR? . In order for the phase transition to occur at the same parameters in 
the limited system as in a bulk one, it is sufficient that the missing potential ener gy is much 
smaller than the kinetic energy. In the case T/c B E , this reduces to R ^> & pC e /Tk B . 
When Tk B <C E , kinetic energy is estimated as h 2 p 2 ^ 3 /m and thus R ^> r (pr 3 l ) 1 / 9 . 

5.3 Results for the phase diagram 

The phase diagram of the model includes a solid at high densities and, at lower densities, 
a gas phase that Bose-condenses at sufficiently low temperatures [OAL + ll]. To locate 
these phase regions, we employed a number of methods, each suitable in a certain area 
of the phase diagram. At zero temperature, the model was treated with the diffusion 
Monte Carlo (DMC), a projector method which provides an exact ground-state energy for 
bosonic systems (see Section 2 of this Thesis). DMC has been used successfully in the past 
to calculate the equations of state and locate quantum phase transitions for a variety of 
systems. Transitions at non-zero temperature were studied with path integral Monte Carlo 
(PIMC) [Gil90, Cha97, Cep95, SCB09], a first principles method which allows to compute 
the averages of quantum operators by summing over the quantum partition function of 
the system. Both DMC and PIMC methods allow to treat systems with several hundred 
particles under periodic boundary conditions, with thermodynamic limits obtained by a 
suitable extrapolation. Additionally, classical limits were established with classical Monte 
Carlo calculations. In two regimes the location of phase transitions could be expressed 
in a semi-analytical form. In the first case, the transition between superfluid and normal 
gas was expressed in terms of the scattering length of the potential by means of a known 
relationship. In the second, the solid-to-gas transition was located at low temperatures 
with the harmonic theory. The results are summarized in the phase diagram shown in 
Fig. 5.3. 

At sufficiently high density, the atoms are expected to form a crystalline solid. Sum- 
ming the potential energy of the perfect lattice structures, we conclude that the preferred 
symmetry is fee. While other structures may be excluded on the energetic grounds, the 
energy of the hep structure is very close to that of the fee. The difference between the per- 
fect crystal energies, E hcp — E {cc = 2 x 10~ 4 (prQ) 2 £ , , is small enough to be comparable to 
or even swamped by the temperature effects in present experiments (for example, in works 
[TFS + 04, HRB + 08, SGH + 10]). The hep phase is anticipated to be metastable with respect 
to the transition to the fee phase. Zero-point motion and temperature effects are expected 
to keep the fee symmetry preferred to hep. If the dressed interaction [LHL12, HNP10] in 
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the form V g - g = l/(£o + r 6 ) between the ground-state atoms is considered (£ stands for a 
blockade radius), the preferable crystal packing can change, as shown in Fig. (5.1). 




-2x10 - 



Figure 5.1: Dependence of the difference of the Madelung energies of ideal fee and hep 
crystals in the units of E Q on blockade radius £o- 

In the rest, we assume the system crystallizes in the fee structure. 
Investigation on the zero-temperature line was done with the DMC method [HJR94, 
BC94]. For importance sampling in the gas phase we used a Jastrow form 

Mr) = LEW [-1/2(V^) 2 ] + exp [-1/2(6/(1 - r 4 ,)) 2 ] } , (5.4) 

Tij < L/2, for a periodic box of size L. The second power in 1/r arises from the cusp condi- 
tion of the scattering problem with the repulsive 1/r 6 potential and is also compatible with 
the presence of long- wavelength phonons [RC67]. The parameter b was variationally opti- 
mized beforehand. The Nosanow-Jastrow wave function was used for importance sampling 
in the solid phase [Nos64, CB08b]. It consists of the product of the above Jastrow term 
and a site-localizing Nosanow term ]li ex P [ — ( r i ~~ h) 2 /^}, where T*j and li denote corre- 
spondingly the coordinates of the atoms and lattice sites, and 7 is the second optimized 
parameter (for the detailed discussion see Section 2.4.2). The breaking of exchange sym- 
metry between particles in the solid affects the energy only negligibly [CACB09]. Within 
the statistical errors of the DMC, results for the energies of the fee and hep lattices are 
indistinguishable and both are lower than the energies derived using bec configuration. 



5.3. Results for the phase diagram 



101 



While the phase transitions are conventionally reported as a function of pressure rather 
than density, density of the Rydberg atoms is more accessible and controllable experimen- 
tally. We therefore choose to express the transition locations in terms of density, even for 
the first-order solidification transition (in this case one needs to specify the coexistence 
region). We find that the equations of state for the fee solid and gas phases cross at the 
transition density 

p c r'l = 3.9 ±0.2, (5.5) 

expressed in the reduced units with the help of Eq. (5.3). The coexistence region of the 
solid and gas phase at zero temperature, determined using the double-tangent Maxwell 
construction, is narrow and is in fact smaller than the above error for the transition den- 
sity (which arises mostly from the extrapolation to the thermodynamic limit; calculations 
were performed with up to 256 particles). The double-tangent Maxwell construction 1 is 
a standard procedure to describe the phase transitions. It relies on equality of chemical 
potentials of two phases at the same pressure in the coexistence region. Alternatively it 
can be presented as a linear dependence of the Helmholtz free energy A on the volume 
V — 1/p in the coexistence region: 

dA x A 2 -A l 

W {Vl) = M (5 - 6) 

and a demand of a constant pressure 

that is equivalent to finding of a common tangent of the equations of two states (at T = 
free energy is equal to the total energy of the system). The Maxwell construction at zero 
temperature is given in Fig. (5.2). The lines cubic polynomial fits to the Monte Carlo data. 
The 

Although the coexistence region is relatively narrow, its width can be obtained with a 
satisfactory accuracy. Let us explain how it is done in our case. First of all we need to 
evaluate the effect of the statistical uncertainty on these results. The Monte Carlo data for 
the equation of state of the liquid and crystalline phases, with the statistic errors for each 
point (it is convenient to work with volume V = 1/p as an argument), is approximated 
with a cubic polynomial. On the first step, the value of energy at each point is moved by 
a random Gaussian shift with the variation, equal to the corresponding statistical error. 
Then, the coefficients of the optimal polynomial approximations for En qu id(V) and E so \ia(V) 
are recalculated, and therefore one finds new positions of the melting and crystallization 
points, the intersection point and the width of the coexistence zone. If this procedure is 
carried out many times, one can obtain the variance of each of the values. In this case the 
error is approximately 1.5%. A separate source of the error in the estimation of the width 



1 For a detailed discussion on the topic see, for instance, [Hua05] 
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Figure 5.2: Study of double-tangent Maxwell construction. Equations of states E(V) (the 
volume V = 1/p, dimensionless) calculated with DMC for the system of 256 particles 
for liquid (red solid line) and solid (green dashed line) phases of the Rydberg system, 
normalized by the intersection energy E . The common tangent is a blue dotted line, blue 
circles stand for two tangent points. The upper and lower plots demonstrate the same data 
with diffenerent ranges of the inverse density. 
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of the coexistence zone is the finite-size dependence of the results, that may be estimated 
by extrapolating the data to N p — > oo. The relative error of this estimation is of the same 
order, as for the phase transition density, given above, that is ~ 5%. The overall result for 
the width of the coexistence region can be evaluated as 



The transition line between solid and gas phases at small temperatures can be de- 
termined with the harmonic theory [AM76], assuming the Lindemann ratio remains un- 
changed on the transition line. The value of the Lindemann parameter at melting may be 
extracted from the DMC calculations of the transition density at zero temperature. The 
resulting low-temperature dependence of the gas-to-solid transition density is given by 



where p c is the transition density at zero temperature, Eq. (5.5), and the constant C = 8.0 
is determined numerically from the dispersion curves of the solid and depends on the 
interactions and geometry of the fee lattice. 

A quantum solid melts at lower temperatures than the classical one due to the zero- 
point motion of the atoms. The classical transition was located in the canonical ensemble 
by Metropolis sampling of the Boltzmann factor. As the potential energy C^/r 6 is exactly 
proportional to the square of the density, the transition temperature for the classical system 
also scales exactly as T oc p 2 . We find that 



As expected, such scaling removes the Planck constant from the classical transition tem- 
perature, which in fact simplifies to J' c classical / CB = 0.22p 2 C*6. 

To fully account for quantum effects, the gas-to-solid transition at T ^ was also 
located with PIMC calculations. We used a decomposition of the action operator that is 
accurate beyond the fourth order [Chi04]. For details of the method and implementa- 
tion, see Ref. [SCB09]. The transition was located by observing melting or solidification 
while working in the canonical ensemble, beginning with configurations of atoms placed 
on a randomly distorted lattice. Used in this way, the calculations determine a range in 
which the transition density is located. PIMC results confirm the validity of the harmonic 
approximation at low temperatures. At higher temperatures the transition density follows 
the classical melting curve (5.10). 

The above results establish the solidification transition of the repulsive van der Waals 
model. Additionally, the dynamic nature of the Rydberg gas raises a possibility for the 
spatial ordering to be induced kinetically, as the combination of decay and strong blockade 
will favor supplanting excitations to be equidistant from their immediate neighbors. We 



Ap = 0.056 ±0.003. 



(5.8) 




(5.9) 




(5.10) 
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modelled such a process and observed that replacement of decaying excitations in local 
blockade regime indeed creates a short-distance order, but not a true long-distance crys- 
talline ordering. These finding are consistent with much more elaborate dynamic models 
of Refs. [RH05, AGHW10, SC10]. 

At low temperature, the gas phase of the model is expected to form a Bose-Einstein 
condensate (BEC). Transition between the BEC and normal gas phases at low densities 
lies slightly above the ideal Bose gas condensation temperature, 



due to the repulsive interaction between particles [PGP08]. The correction is governed 
by the scattering length of the potential a s , which can be found to be equal to a s = 
2 r (3/4) /r (1/4) 7*0 = 0.676 . . . r . The transition temperature is then given by Tbec = 
t bec (l + ca s p 1/3 ) , where c is a positive constant of the order of unity (for details, see 
Ref. [PGP08] and references therein). In the present case this expression is only valid at 
very low densities (one needs to satisfy at least pr^ < 5 x 1CT 2 to make the description in 
terms of the zero- momentum scattering length meaningful), where the magnitude of the 
correction is not significant. 

At higher densities the BEC-to-normal gas transition is no longer universal and depends 
on the form of the potential. We determined the location of this second-order transition 
with the PIMC method by calculating the superfluid transition from the winding number 
estimator [PC87b]. The PIMC calculations show that at higher densities the interactions 
deplete the condensate and the transition temperature is lower than for the ideal Bose 
gas. Combining the PIMC results, the region in which the triple point is located was 
determined as 4.5 < T/(Eo/k&) < 6.5 and 4 < pr^ < 5, which we consider sufficiently 
narrow for practical considerations. 

5.4 Comparison with experimental conditions 

Because the interaction constant C 6 enters the reduced units [Eq. (5.3)], the effective 
temperature and density can be varied over many orders of magnitude. Most of the present 
experiments are deeply in the "classical" region of the phase diagram (Fig. 5.3). As an 
example, we consider the conditions of the experiments presented in Ref. [HRB + 08]. For 
the excitation with 170 ns laser pulses, the system parameters at 4 pK are T/(£'o//cb) ~ 
33 x 10 5 and pr$ ps 1.9 x 10 3 , which in fact correspond to the gas phase of the equilibrium 
phase diagram. For 320 ns excitation pulses and T = 1 pK, T/(E /k-B) ps 8.2 x 10 5 
and prQ ps 7.4 x 10 3 , well below the gas-to-solid transition. Therefore, the achievable 
temperature and density are already in the range suitable for investigating the equilibrium 
phase diagram. Increasing the excitation number increases the interaction constant Cq 
and moves the system deeper into the classical regime where the gas and solid phases are 
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Figure 5.3: Phase diagram of the repulsive C^/r 6 interaction, scaled to units given in 
Eq. (5.3). Location of the gas-to-solid transition at zero temperature, determined with 
DMC, is shown with an open blue square on the T = axis. Dashed blue line shows 
gas-to-solid transition as found with harmonic theory [Eq. (5.9)]. Solid blue line shows 
classical gas-to-solid transition [Eq. (5.10)]. Solid blue squares show location of the gas-to- 
solid transition as determined with PIMC. Red short-dashed line marks the Bose-Einstein 
condensation of the ideal gas [Eq. (5.11)]. Filled red bullets show Bose-Einstein conden- 
sation temperatures found with PIMC. 
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separated by the simple condition of Eq. (5.10). The quantum regime of the phase diagram 
may be accessed by decreasing the excitation numbers or increasing the Forster defect 5. 

Whether Rydberg atoms in actual experiments will reach or even approach an equi- 
librium phase depends on their lifetime, the experiment's duration and availability of a 
relaxation mechanisms. Because of the short lifetimes of the Rydberg states, most current 
experiments are performed on such short timescales as to make the thermal motion negli- 
gible. It is therefore said that the experiments are performed with Rydberg excitations of 
a frozen gas. If the experiments are extended closer to the currently achievable lifetimes of 
the Rydberg states, which can be as large as 100 /is [vDKH + 08, LWK + 09], some degree 
of thermal equilibration will already be achieved. Besides the thermal motion there are, 
however, at least two other kinds of motion that may need to be considered. The first one 
is the motion of the excited atoms due to the strong forces between them. The character- 
istic timescale associated with such a motion is the time that it takes for a Rydberg atom 
to travel the mean distance between Rydberg atoms. Given the mean distance £ ~ p -1 / 3 
and the imbalance force of the order of C^/^ 6 , this time is given by 



ballistic ~ y ' ^' 12 ^ 

For a small fixed number of Rydberg excitations time (5.12) decreases rapidly with the 
excitation number n as n -11 / 2 ; as the local blockade is reached, p oc Cq 1 ^ 2 , and ^ballistic 
instead grows as n 11//6 . For example, Rydberg systems created by the 1.970 /is pulses 
from \\iK gas in the experiment of Ref. [HRB + 08] have tbaiiistic ~ 12 /is. For the setup 
of Ref. [LWK + 09], ^ballistic ~ 60 /is while the clouds could be successfully studied for as 
long as 20 ps. Collisional ionization and heating could potentially hamper such relax- 
ation [ARLGW07]. 



5.5 Excitation spectrum of a classical crystal 

In the previous parts of this Section the properties of a quantum system comprised of 
Rydberg atoms at zero and finite temperature were investigated, and the quantum phase 
diagram was presented. 

In the present subsection, we consider a classical system, governed by the van der Waals 
interaction 1/r 6 and look for its excitation spectrum in the low temperature regime (or 
equivalently in the regime of very strong interactions) . In accordance with the comparison 
of classical potential energies of different kinds of crystal lattice, for T = we chose the 
face-centered cubic (fee) lattice as preferable, with the Madelung energy of the hexagonal 
close pack (hep) being only slightly higher. In practice the fee packing is also advantageous, 
since its elementary cell is "cubic", and its first Brillouin zone has central, plane and axis 
symmetry. 1 The first Brillouin zone of fee formation is known to be a truncated octahedron 

lr The first Brillouin zone is a primitive cell, not reducible to a smaller one by translations of the lattice 
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with the critical points K,L,U,W,X (see Fig. 5.4). 




Figure 5.4: Brillouin zone of fee crystal packing. Figure is taken from Wikipedia webpage 
http:/ / en. wikipedia. org/wiki/Brillouin_zone. 



The Hamiltonian of the system takes the form 



N 



a 



i=l i>j I ' » J I 

We consider the case of low kinetic energy, that is when the interaction strength constant Cq 
is large enough, and the system may be treated in a perturbative manner within a harmonic 
approximation. The positions of the crystal nodes can be fixed at In order to find the 
excitation spectrum of this system we can follow the procedure of Refs. [BM77, MPW07]. 
First of all, we give to the particles small arbitrary displacements Ui from the nodes. The 
potential energy is then represented by 

E P = ^J2] " T ( 5 - 14 ) 

^ ijLj \ r ij + U ij\ 

with the notations = r.i — rj, = — Uj and k = 6. 



in reciprocal space 
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Then, we perform an expansion of the potential energy E p in powers of the displace- 
ments uf, where a labels the Cartesian coordinates a = x,y,z and U stands for the 
Madelung energy of a chosen crystal packing. 

e p = u + \t.T. + - ( 5 - 15 ) 

ia j/3 

The first-order term here vanishes due to the fact that {Ri} is a minimum-energy config- 
uration. The Hessian matrix $ Q( g(rjj) is obtained through a simple double differentiation: 

C k (k(k + 2)r°4 kS aP , 

if l^], (5.16) 

(5.17) 

The Fourier transformed dynamical matrix C a /3(q) is related to $ Q/ g (ry) by the expres- 
sion 

C a p{q) = — ®*p( r ij) exp[-iqnj] (5.18) 

where m stands for the mass of a particle in the kinetic term of the Hamiltonian. In order 
to simplify the treatment and get rid of the conditional definitions for $ QJ g, we introduce 
the matrix S a p(q) such that 

^( k(k + 2)rg-rg k5ap\ . 
SaM) = L fc+4 - "TO expf-igry] (5.19) 





then it can be seen that 



Capiq) = -7^(Safi(q) - SapiO)) (5.20) 

The dynamical matrix can be used to obtain the normal mode frequencies by solving the 
eigenvalue problem 

C a p(q)e(q) = co 2 e(q). (5.21) 

Since C a p is a real symmetric 3x3 matrix, there exists a complete and orthonormal set of 
three eigenvectors, written in which the Hamiltonian is diagonal. The wave vectors q are 
taken from the irreducible Brillouin zone. 

For our case of the chosen model potential between Rydberg atoms k — 6, one has to 
evaluate the following sums 
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The Hamiltonian is rewritten in dimensionless form as 

s i=l i<j I 1 j\ 

where distance and energy units are 

ro = l/p 1/3 (5.25) 

E = C 6 /r 6 , (5.26) 

and the following form for the dimensionless parameter r s 

l h 2 h 2 



r s mrQ-Eo mCgp 4 / 3 ^ 

which characterizes the ratio between the kinetic and potential energies. It can be noticed 
that the parameter r s plays the role of the mass m in Eq. (5.18). The units of frequency 
of the normal modes are taken as 

1 h 

W ° " VF= " (mC 6 )VV/ 3 • (5 ' 28) 
In reduced units the problem can be rewritten as 

u 2 (q)u = <5> ctP (q)u (5.29) 

with 

The kernel $> a p(q) can be calculated by means of the Ewald summation technique. 

The excitation spectrum in this units is shown in Fig. 5.5. In our description the three 
distinct branches of the spectrum can be degenerate (the second and the third coincide on 
the path GX) . For small values of momentum q the excitations are in phononic (linear) 
regime modes with the sound velocity c = u}/\q\, which is dependent on a direction in a 
crystal. For GX direction in the crystal two sound velocities (of two coinciding transverse 
modes and one longitudinal mode) are equal to 

ci = 19.41 
c 2i3 = 10.04 

in the units of u)o r o- It can be noticed that the longitudinal mode is also phononic unlike 
the case of the three-dimensional Wigner crystal, where a typical long-wavelength behavior 
of the frequency is 0Ji(q) ~ y/q. The spectrum is periodic as one reaches a border of the first 
Brillouin zone, the values of the frequencies are also continuous, as seen in Fig. 5.5, when a 
momentum vector approaches the point X from different directions. At finite temperature 
the spectrum is expected to be smeared. 
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q[27i/r ] 

Figure 5.5: Excitation spectrum of the van der Waals 1/r 6 interaction. The wave vector q 
follows the path G,X,K,L,U,W,X in the first Brillouin zone with the metrics of the vector 
preserved. The units for frequency and wave vector are loq and 27r/r respectively, see 
Eqs. (5.25) and (5.28). 

5.6 Discussion and conclusions 

A very interesting and open question is the possibility of supersolidity in Rydberg systems. 
Ground state atoms dressed in Rydberg states exhibit weak van der Waals interactions at 
large distances, as described in [HNP10]. The question of supersolidity of such atoms was 
addressed in Refs. [CJB + 10, HNP10]. Here, we consider an alternative scenario in which 
the gas is additionally allowed to have a lattice of Rydberg excitations. Such a lattice would 
in turn impose weak but long-ranged spatial correlations onto the ground state atoms. At 
the same time, the ground-state atoms may be Bose condensed [HRB+08]. However, it 
may be impossible to identify which of the atoms was excited within a certain proximity, as 
was demonstrated, for example, by the superatom analysis of the experimental results in 
Refs. [HRB+07, TFS+04]. However, motion will lead to dephasing of this state [HLW+11]. 
If the atoms are indeed prepared in such a mixed state, combining the ground \g) and 
excited |e) states as \g) + a\e), N" 1 < a C 1, then both the lattice-forming and the 
BEC components are indistinguishable and may be said to be formed by the same atoms. 
Therefore, such a system would consist of particles which would simultaneously break 
translational symmetry and possess off-diagonal long-range order, which is a realization of 
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supersolid. While our model does not include the light field, the conditions for the phases 
of both excited and ground state atoms may be immediately extracted from Fig. 5.3, just 
with different reduced units for the two species. 

In conclusion, it is possible to parameterize a model with isotropic van der Waals in- 
teractions into a universal phase diagram. We have characterized the phase diagram of 
Rydberg atoms by considering a model of bosons with repulsive van der Waals interac- 
tion, and determined solidification and Bose-Einstein condensation conditions. Relaxation 
mechanisms other than thermal motion should be considered if one considers Rydberg sys- 
tems on timescales of several tenth of microseconds. We have also studied the excitation 
spectrum within the approximation of a classical harmonic crystal. Finally, it is worth 
mentioning that interactions between Rydberg excitations open a possibility of new super- 
solid scenarios. 



Chapter 6 

Para-hydrogen at low temperature 



6.1 Introduction 

Superfluidity and Bose-Einstein condensation (BEC) have been stunningly shown in me- 
tastable dilute alkali gases, magnetically confined at ultralow temperatures. [PS03] The 
extreme diluteness of these gases allows for the achievement of BEC with an almost full 
occupation of the zero-momentum state that has been possible to observe and measure 
quite easily. This contrasts with the difficulties encountered in the measure of the conden- 
sate fraction in liquid 4 He, which amounts only 8% at the equilibrium density. [GDA + 11] 
However, liquid 4 He is a stable superfluid below the lambda transition T\ = 2.17 K and 
therefore a system more easily accessible. Before the blowup produced in the field of quan- 
tum fluids by the first experimental realization of BEC gases, liquid helium was the only 
paradigm of a superfluid. From long time ago, there has been great interest in the search 
of superfluid condensed phases other than liquid helium. Spin-polarized atomic deuterium 
and tritium are predicted to be fermionic and bosonic liquids, respectively, in the limit 
of zero temperature. [PC87a, BMVB09] However, its experimental study has proven to 
be very elusive due to its high recombination rate, and only the case of atomic hydrogen, 
whose ground state is a gas, has been experimentally driven to its BEC state. [FKW + 98] 
The next candidate for superfluidity is molecular hydrogen, which has been studied for a 
long time. [Sil80] This seems a priori an optimal system due to its very light mass but it 
crystallizes at relatively high temperature as a consequence of the intensity of its inter- 
molecular attraction, without exhibiting any superfluid transition in the liquid phase. In 
the present work, we study the properties of metastable liquid or glass molecular hydrogen 
at very low temperatures using quantum Monte Carlo methods. 

In 1972, Ginzburg and Sobyanin [GS72] proposed that any Bose liquid should be su- 
perfluid below a certain temperature T\, unless it solidifies at temperature Tf higher than 
T\. To give a first estimation of T\, they used the ideal Bose gas theory, obtaining 

where m is the atomic mass, g is the spin degeneracy, ks is the Boltzmann constant and 
p is the density of the system. Ginzburg and Sobyanin proposed molecular para-hydrogen 
(p-H 2 ) as a plausible candidate for superfluidity: being a spinless boson (g = 1) with a 
small mass, p-B.2 should undergo a superfluid transition at a relatively high temperature 
(according to Eq. (6.1), T x — 6 K). 
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The estimation of T\, given by Eq. (6.1), is clearly inaccurate in the case of dense 
liquids because it cannot account for the observed dependence of T\ with the density. In 
fact, T\ slightly decreases in liquid 4 He when p increases, a manifestly opposite behavior 
to the increase with p 2 / 3 given by the ideal gas formula (6.1). In order to provide a more 
reasonable estimation of T\, Apenko [Ape99] proposed a phenomenological prescription for 
the superfluid transition, similar to the Lindemann criterion for classical crystal melting. 
In this way, he was able to take into account quantum decoherence effects due to the strong 
interatomic potential and to relate the critical temperature for superfluidity with the mean 
kinetic energy per particle above the transition. For p-H 2 , he concluded that T\ should 
vary between 1.1 K and 2.1 K, depending on the density of the system. 

Superfluid p-H.2 is not observed in a stable form because it crystallizes at temperature 
Tf = 13.8 K, which is significantly higher than the expected T\. Several studies about 
crystal nucleation in p-H 2 have been performed in order to understand if the liquid can 
enter a supercooled phase, i.e., a metastable phase in which the liquid is cooled below its 
freezing temperature without forming a crystal. Maris et al. [MSH83] calculated the rate 
r(T) of homogeneous nucleation of the solid phase from the liquid as a function of the 
temperature T, showing a maximum of T around T = 7 K and a rapid decrease at lower 
temperature. This suggests that, if it would be possible to supercool the liquid through 
the range where T is large, one might be able to reach a low-temperature region where 
the liquid is essentially stable. However, recent experiments have indicated that, even at 
T~9K, the rate of crystal growth is so high that the liquid phase freezes quickly into a 
metastable polymorph crystal. [KFT + 11] 

Even though several supercooling techniques have been proposed to create a metastable 
liquid phase in bulk p-H.2, [MSW87, VMOO, GFP + 06] none of them has proven so far to 
be successful and no direct evidence of superfluidity has been detected. However, there 
are evidences of superfluidity in several spectroscopic studies of small doped p-H 2 clusters. 
In 2000, Grebenev et al. [GSTV00] analyzed the rotational spectra of a linear carbonyl 
sulfide (OCS) molecule surrounded by 14 to 16 p-H 2 molecules absorbed in a larger helium 
droplet, which fixes the temperature of the cluster. When P-H2 is immersed in a 4 He 
droplet (T = 0.38 K), the measured spectra shows a peak indicating the excitation of 
angular momentum around the OCS axis. On the other hand, if the small p-H 2 cluster is 
put inside a colder 4 He- 3 He droplet (T = 0.15 K), the peak disappears: the OCS molecule 
is then able to rotate freely inside the hydrogen cluster, pointing to the superfluidity of 
the surrounding P-H2 molecules. These results have been confirmed in a later experiment 
on small P-H2 clusters doped with carbon dioxide. [LLRRM10] From a precise analysis of 
the rotational spectra, it has been possible to measure the effective momentum of inertia 
of these small systems, and thus of their superfluid fraction p s , providing a clear evidence 
of superfluidity in clusters made up of N < 18 p-H 2 molecules. These clusters are too 
small for extracting reliable predictions of a metastable liquid phase and larger clusters 
would be desirable. To this end, Kuyanov-Prozument and Vilesov [KPV08] have been 
able to stabilize liquid clusters with an average size of N rs 10 4 p-H 2 molecules down to 
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temperature T = 2 K, but they do not see any evidence of superfluidity. Other attempts of 
producing liquid p-H 2 well below Tf (T — 1.3 K) are based on the generation of continuous 
hydrogen filaments of macroscopic dimensions. [GFP + 06] 

The search for a superfluid p-H 2 phase has been intense also from the theoretical point of 
view. The rather simple radial form of the £>-H 2 -j>-H 2 interaction and the microscopic accu- 
racy achieved by quantum Monte Carlo methods have stimulated a long-standing effort for 
devising possible scenarios where supercooled p-B. 2 could be studied. In practically all the 
cases, the search is focused on systems of reduced dimensionality or in finite systems. PIMC 
simulations of p-H.2 films adsorbed on a surface with impurities observed superfluidity for 
some arrangements of these impurities, [GC97] but these results were posteriorly questioned 
by other PIMC studies. [Bon05] In a one- dimensional channel, like the one provided exper- 
imentally by narrow carbon nanotubes, it has been predicted a stable liquid phase in the 
limit of zero temperature. [GBCOO] The largest number of theoretical works are devoted to 
the study of small clusters, both pure [SCK91, MB06, MB07, KSCT07, MB08, SB11, GN08] 
and doped with impurities. [KW02, PZKW05, KW05] All these simulations show that p- 
H 2 becomes superfluid below a certain temperature T = 1-2 K and that the superfluid 
fraction depends on the number of molecules of the cluster. When the cluster becomes 
larger than a certain molecular number (N > 18-25), solid-like structures are observed and 
the superfluidity vanishes. 

We address the calculation of the equation of state of the metastable liquid p-H. 2 phase 
in the limit of zero temperature using the diffusion Monte Carlo (DMC) method. The 
simulation of the liquid phase in this limit is easier than at finite temperature and therefore 
it is able to provide accurate information on its main energetic and structure properties. 

The rest of the chapter is organized as follows. In Sec. 6.2, we introduce the quantum 
Monte Carlo methods used in the study and report specific details on how the simulations 
are carried out. Sec. 6.3 contains the results of the equation of state, structure properties, 
and condensate fraction of metastable liquid p-H 2 at zero temperature, and finally the 
main conclusions of the present work are discussed in Sec. 6.4. 

6.2 Methodology and construction of trial wave func- 
tions 

The H2 molecule, which is composed of two hydrogen atoms linked by a covalent bond, 
is spherically symmetric in the para-hydrogen state (total angular momentum zero). The 
energy scale involved in electronic excitations (~ 10 5 K) is orders of magnitude larger than 
the intermolecular one (~ 10 1 K), thus to model the p-H 2 -p-H 2 interaction by means of a 
radial pair-potential and to consider the molecules as point-like turns out to be justified 
upon the condition of low or moderate pressures. In this work, we have chosen the well- 
known and commonly used semiempirical Silvera- Goldman pair potential. [SG78] This 
potential has proved to be accurate at low temperature and at the pressure regimes in 
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which we are interested. 

The study in the limit of zero temperature has been performed with the DMC method. 
DMC is a first-principles method which can access exactly the ground state of bosonic 
systems. It is a form of Green's Function Monte Carlo which samples the projection of the 
ground state from the initial configuration with the operator exp [— (H — E )t}. Here, % 
is the Hamiltonian 

fc2 N N 

^ = -5rE V '+ E V( rij ), (6.2) 

z " t i=l l=i<j 

E is a norm-preserving adjustable constant and r is the imaginary time. The simulation 
is performed by advancing in r via a combination of diffusion, drift and branching steps on 
walkers -R (sets of 3N coordinates) representing the wavefunction of the system. [BCN94] 
The imaginary-time evolution of the walkers is "guided" during the drift stage by a guiding 
wavefunction (fro, which is usually a good guess for the wavefunction of the system. This 
function contains basic ingredients of the system as its symmetry, phase and expected 
behaviors at short and long distances according to its Hamiltonian. Technically, (frc al- 
lows importance sampling and thus reduces the variance of the ground-state estimations. 
It is straightforward to show that for the Hamiltonian % and any operator commuting 
with it, the expectation value is computed exactly within statistical error. Other diago- 
nal operators which do not fulfill this condition require of a special treatment, known as 
pure estimation, [CB95b] which leads also for this case to unbiased results (for details see 
Section 2.5 of this Thesis). 

The phase of the system is imposed within the typical imaginary-time length by the 
guiding wave function. This property of the DMC method is here a key point if we are 
pursuing a prospection on the properties of the metastable liquid p-B. 2 phase. Then, for 
the liquid phase <Pg is taken in a Jastrow form 

N 

<P G (R) = I] f( rij ) , (6.3) 
l=i<j 

with a two-body correlation function [Rea79] 

2" ~ 



f{r) = exp 




r — A 
A 



(6.4) 



In order to compare the results obtained for the liquid phase with the ones corresponding 
to the stable hep solid we have carried out some simulations with a guiding wave function 
of Nosanow- Jastrow type 

N N 

th(R)= II /(%) lidiru) , (6.5) 

l=i<j i=l 

the set {ri} being the lattice points of a perfect hep lattice. Optimal values for the 
parameters entering Eq. (6.4) are b = 3.68A, L = 0.2, A = 5.24A, and A = 0.89A for the 
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liquid phase, and b = 3.45A, L = 0.2, A = 5.49A, and A = 2.8lA for the solid one. The 
Nosanow term is chosen in Gaussian form, g(r) = exp(— 7r 2 ). The density dependence of 
the parameters in the Jastrow term is small, and neglected in practice when used in DMC, 
whereas the Nosanow term parameter 7 is optimized for the whole range of densities. We 
have used 256 and 180 particles per simulation box for the liquid and hep solid phases, 
respectively. The number of walkers and time-step have been adjusted to reduce any bias 
coming from them to the level of the statistical noise. 

6.3 Results at zero temperature 

We have calculated the main properties of the metastable liquid and stable hep solid 
phases of p-H.2- Our main goal has been to know the properties of a hypothetical bulk 
liquid phase and compare them with the ones of the stable solid. In order to achieve reliable 
estimations of liquid p-H 2 it is crucial to work with a guiding wave function of liquid type, 
as we have discussed in the preceding Section. Within the typical imaginary-time length of 
our simulations we have not seen the formation of any crystal structure, i.e., no signatures 
of Bragg peaks in the structure function S(k) have been registered so far. 

In Fig. 6.1, we plot the DMC energies per particle of metastable liquid p-H 2 as a 
function of the density. For comparison, we also report the results obtained for the hep 
crystal phase. Our hep energies are in close agreement with the ones reported in Ref. [OP06] 
using the same Silver a- Goldman potential. In the figure, we also show the experimental 
estimation at T = K from Ref. [Sch70], E/N = —89.9 K, that lies a bit below of our 
results. This is again in agreement with previous DMC results [OP06] which show that 
the experimental energy is, in absolute value, underestimated and overestimated by the 
Silvera-Goldman and Buck potential, [BHK + 83] respectively. Our results for both phases 
are well reproduced by the polynomial law 



(E/N)o and po being the equilibrium energy per particle and equilibrium density, respec- 
tively. These equations of state are shown in Fig. 6.1 with lines. The optimal parameters of 
the fits are: p = 0.026041(20)A" 3 , (E/N) = -86.990(37) K, A = 232(2) K, B = 156(11) 
K for the solid, and p = 0.023386(40)A" 3 , (E/N) = -76.465(51) K, A = 188(1) K, 
B = 131(10) K for the liquid. As expected, our DMC results shows that the solid phase is 
the stable one with a difference in energy per particle at the respective equilibrium points 
of ~ 10 K, the equilibrium density of the liquid being ~ 10 % smaller than the solid one. 
The same trend was observed in a DMC simulation of two-dimensional p-H.2, but there 
the differences were significantly smaller. [CB08a] It is worth noticing that about one half 
of the energy difference in the bulk systems comes from the decrease of the kinetic energy 




(6.6) 
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Figure 6.1: DMC energies per particle as a function of the density. Squares and circles 
correspond to the liquid and solid phases, respectively. Solid and dashed lines are the 
polynomial fits to the DMC energies for the liquid and solid, respectively. The diamond is 
the experimental energy of hep molecular hydrogen from Ref. [Sch70] 
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per particle going from the liquid to the solid: at density p = 0.03 A 3 , it amounts 93.3(1) 
and 89.1(1) K for the liquid and solid, respectively. 

From the equations of state (6.6), it is easy to know the pressure of the system at any 
density using the relation P(p) = p 2 [d(E/N)/dp\. The results obtained for metastable 
liquid and stable solid phases are shown in Fig. 6.2. As one can see, at a given density 
the pressure of the liquid is larger than the one of the solid mainly due to the different 
location of the equilibrium densities (P = 0). The results for the solid are compared 
with experimental data from Ref. [DdWS79] showing a good agreement especially for not 
very large pressures. The density at which the function P(p) has a zero slope defines the 
spinodal point; beyond this limit the system is no more thermodynamically stable as a 
homogeneous phase. At this point, the speed of sound c(p) = [m~ l (dP / dp)] 1 / 2 becomes 
zero. Results for c(p) are shown for both phases in Fig. 6.3. The speed of sound decreases 
when the density is reduced and drops to zero at the spinodal point: (p c = 0.0176(1) A , 
P c = -12.6(5) MPa) and (p c = 0.0195(1) A^ 3 , P c = -17.5(5) MPa) for liquid and solid, 
respectively. 

DMC produces also accurate results for the structure of the bulk system. In Fig. 6.4, we 
show results for the two-body radial distribution function g(r) of the liquid p-B.2 phase for 
a set of densities. This function is proportional to the probability of finding two molecules 
separated by a distance r. Increasing the density, the main peak becomes higher and moves 
to shorter interparticle distances; at least three peaks are observed. All these features point 
to the picture of a very dense quantum liquid, with much more structure than in stable 
liquid 4 He. In the same Fig. 6.4, we show results for the static structure factor S(k), 
related to g(r) by a Fourier transform. As one can see, the main peak increases quite fast 
with the density suggesting a highly structured metastable liquid. Nevertheless, we have 
not observed within the scale of the simulations the emergence of any Bragg peak which 
would point to formation of crystallites in the simulation box. In Fig. 6.5, we illustrate the 
comparison between S(k) for the liquid and solid systems at a density p = 0.0245 rA -3 , 
close to the equilibrium density of the liquid. The difference is the one expected between a 
liquid and a solid: oscillating function towards one at large k for the liquid and a sequence 
of Bragg peaks, corresponding to the hep lattice, for the solid. 

One of the most relevant properties of a superfluid is the mean occupation of the zero- 
momentum state, i.e., the condensate fraction no- As it is well known, no can be obtained 
from the asymptotic behavior of the one-body density matrix p(r), 

n = lim p(r) , (6.7) 

r->oo 

with p(r) being obtained as the expectation value of the operator 

/ $(ri,...,r i + r,...,r JV ) \ 

\ *(r„...,r w) / ' (6 - 8) 

DMC results for the condensate fraction of liquid p-B.2 as a function of the density, obtained 
using the extrapolated estimator (there are no reliable pure estimators for non-diagonal 
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operators), are shown in Fig. 6.6. The decrease of uq with the density is well described 
by an exponential decay (line in the figure). The strong interactions induced by the deep 
attractive potential well produce a big depletion of the condensate state. At the equilibrium 
density, our estimation for the condensate fraction is n = 0.0037(7). This value is more 
than one order of magnitude smaller than the measured condensate fraction [GDA + 11] of 
liquid 4 He at equilibrium (0.08). 

6.4 Conclusions 

A possible observation of superfluidity in bulk or clustered systems of para-hydrogen atoms 
p-H 2 was a subject of an extensive study in the last decades ( [GS72], [SCK91]) due to a va- 
riety of advantageous properties of such atoms. Nevertheless the experimental observation 
has not been performed yet, mostly because of a relatively high solidification tempera- 
ture of p-H.2, that can bring a perspective superfluid transition deeply into the metastable 
phase. 

Our Group performed a multifold study of the system of para-hydrogen atoms p-H 2 at 
low temperatures deeply below the point of crystallization by means of Quantum Monte 
Carlo methods. The zero-temperature simulation was performed in order to investigate the 
properties of a metastable liquid phase and to find the fraction of the Bose-Einstein conden- 
sate in the relevant range of densities. The methods of choice for the zero-temperature sim- 
ulations of the p-H 2 system were the Variational Monte Carlo and the diffusion Monte Carlo 
techniques. The latter is an efficient and versatile instrument to calculate the quantum 
properties of the system, including the non-local ones. The results of the zero-temperature 
simulations suggest that the metastable liquid para-hydrogen is a strongly correlated liq- 
uid, that again might be a sign of instability of this hypothetical system. The calculation 
of the Bose-Einstein condensate show that the condensate fraction is substantially lower 
than in the liquid helium 4 He. 
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Figure 6.2: Pressure of the liquid (solid line) and solid (dashed line) p-H.2 phases as a 
function of the density. Experimental points for the solid phase [DdWS79] are show as 
solid circles. 
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Figure 6.3: Speed of sound of the liquid (solid line) and solid (dashed line) p-H 2 phases as 
a function of the density. 
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Figure 6.4: Top panel: Two-body radial distribution function of the liquid p-H.2 phase at 
different densities: solid, long-dashed, short-dashed, and dotted lines stand for densities 
p = 0.0180, 0.0245, 0.0300, and 0.0340 A , respectively. Bottom panel: Static structure 
factor of the liquid phase. Same densities and notation than in the top panel. 
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Figure 6.5: Static structure function of liquid and solid p-H 2 at density p = 0.0245 A 3 . 
The result for the liquid S^k) (left scale) is shown with a solid line; the one for the hep 
solid Ss{k) (right scale) with a dashed line. 



6.4. Conclusions 
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Figure 6.6: Condensate fraction of metastable liquid p-B. 2 as a function of the density. The 
points are the DMC results and the line is an exponential fit to them. 



Chapter 7 

Conclusions and list of publications 



Here we summarize the results, presented in the Thesis. 

In Chapter 3, devoted to the generalization of the Ewald summation technique, we 
derived explicit expressions for the Ewald sums in three-dimensional systems, governed by 
a generic l/|r*| fc power-law interaction, with periodic boundary conditions applied. We also 
extended the derivation to the cases of two-dimensional and one-dimensional geometry. 
The importance of these generalizations can be seen, as they apply to many physically 
important interaction potentials as the dipole-dipole interaction, van der Waals interaction 
etc. In this Thesis we give the functional forms for the terms of the Ewald sums (both 
in momentum and coordinate space), ready for implementation in actual calculations. 
The derivation and the functional form of the results differs in the cases of short-ranged 
(k > D), long-ranged {k < D) and "marginal" (k = D) forces, where D is the system 
dimensionality. The cases of long-range forces require separate calculations because of the 
divergence of the potential energy. This unphysical divergence disappears when demanding 
a charge-neutrality of the system, which can be restored within the "jellium" model. The 
resulting expressions in the jellium model are explicitly given. It is argued that in the case 
of some short-range potentials the Ewald method can be advantageous with respect to a 
direct summation due to a faster convergence rate, typically, of a Gaussian kind versus 
a certain power-law rate. We also give a discussion of the convergence properties of a 
quasi-neutral Coulomb system. 

The results are first presented for the case of a 3D system in a cubic simulation box 
in order to explain the general mathematical procedure, which for the specific case of the 
Coulomb potential recovers well-known results [AT89]. Later on, the same mathematical 
technique is applied to 2D and ID geometries. For the one- dimensional case the initial 
sum for the potential energy is explicitly evaluated in Eq. (3.96), nonetheless the Ewald 
summation is developed as well for this case and may be used as a mathematical equality. 
The more efficient representations of the reciprocal space sums, which have a complexity 
O(N) instead of 0(N 2 ) with the number of particles N, are presented for 3D and 2D sys- 
tems. The explicit expressions for the terms of the Ewald sums are given in a tabular form 
for physically relevant potentials with small integer power indices k, as dipole-dipole inter- 
action potential, Lennard- Jones potential and others in both three- and two-dimensional 
geometries (refer to Tables 3.1 and 3.2). 

When the simulation box cannot be chosen cubic, for example in a modelling of a three- 
dimensional hep crystal structure, the Ewald method can also be applied after certain 
modifications. Formally, it consists in the choice of an appropriate rectangular simulation 
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box and a substitution of the vector n by n r = (n x L x + n y L y + n z L z )/Lo and = 
{m x /L x + rriy/Ly + m z /L z )L in the real and momentum space sums, respectively [see 
Eqs. (3.102) and (3.103)]. 

The optimization of the involved parameters, that is the integration parameter a and 
the cut-off lengths in coordinate and momentum space, is a necessary operation in order 
to improve the convergence rates and to avoid excessive calculations. The main idea of the 
optimization, that we propose in the this Thesis, is to perform a benchmark calculation, 
minimizing the variance in the energy. We show how this optimization technique works 
on the example of a two-dimensional gas of dipoles, aligned perpendicularly to the plane 
of motion. This proposed optimization technique is mostly efficient in simulations of gases 
and crystals. In spite of being very simple, it allows to find rather quickly adequate 
parameter ranges. The analytical estimations of the optimal parameters are given as well 
and are proven to be consistent with the results, obtained in an explicit optimization. 
A more sophisticated method to optimize the calculation parameters, taking advantage 
of the O(N) representation of the Fourier transform sum, is also presented with explicit 
estimations of the parameters for a typical system simulated by Quantum Monte Carlo 
methods. 

Another problem studied is the phase diagram of Yukawa systems. The Coulomb 
potential is probably one of the most basic and widespread interactions in the nature as 
it describes forces between charges. In a charge neutral system, the presence of a second 
component introduces screening. The Yukawa potential as well appears in the problem 
of mass-imbalanced two-component Fermi gas. It is known that in the case of a large 
mass imbalance, the light fermions introduce an effective Yukawa potential between the 
heavy fermions. It was predicted in [PAP + 07] that this effective interaction might lead 
to crystallization in two-dimensional systems, although no estimations were done in three- 
dimensional case. The study of Ceperley et al. [CCK76] partially addresses the problem of 
the zero-temperature phase diagram. In their work the transition line was drawn by using 
approximate Lindemann melting criterion. For the first time a fully quantum mechanical 
calculation is done in order to find the zero-temperature phase diagram of Yukawa particles. 
The diffusion Monte Carlo method is expected to give a ground state energy exactly, so 
the phase diagram should be exact. The melting and solidification points are found by 
using the double-tangent Maxwell construction. The Lindemann criterion indeed applies 
in a large part of the phase diagram apart from the high density regime where notable 
differences compared to the prediction of Ceperley et al. are visible. 

We discovered that an improvement to the efficiency of the diffusion Monte Carlo 
method can be obtained by using the hypernetted chain (HNC) method (based on a solution 
of the corresponding Euler-Lagrange hypernetted chain equations) in a construction of 
Jastrow guiding wave functions. The HNC method allows to find a very good numerical 
approximation for the two-body Jastrow terms provides a result, which is already optimized 
in a many-body sense. The HNC solutions are also advantageous, since they do not require 
cumbersome multiparametric optimizations, and are quite precise and very fast. 
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The phase diagram of Yukawa particles has a peculiar shape. That is, a double tran- 
sition is possible, when for a fixed interaction strength, a change in the density might 
cause reentrant gas-solid-gas transition. This unusual behavior is caused by a competi- 
tion between long-range Coulomb part of the Yukawa potential (melts at large densities) 
and exponential screening (melts at small densities). To the best of our knowledge this 
is the first time that the high-density melting was observed in this system in a quantum 
simulation. 

Based on our calculations we find that the mass ratios in achievable fermionic mixtures 
of pure elements are too small to undergo a transition to a crystal phase regardless of the 
density of the bulk. Nonetheless, larger effective mass ratios can be reached if the heavier 
component is confined by an optical lattice, that enhances strongly the strength of the 
lattice. This alternative procedure of obtaining large effective mass ratios is discussed in 
the Thesis, and based on the available data on the current and perspective experimental 
set-ups the principal feasibility to produce experimental conditions for observation of the 
crystallization is argued. The heights of optical lattices are typically tunable by the in- 
tensity of laser beams, so that the correctness of our theoretical predictions for the phase 
transition can be checked. 

In order to localize more accurately the phase transition point in the high density 
regime we applied the Ewald summation technique, that yields the potential energy of the 
simulation cell, replicated infinitely in the space. This method proved to be efficient in 
enhancing the rate of the convergence of the DMC energy with the number of particles 
in spite of introducing additional calculation costs. The specific Ewald summation for the 
case of a Yukawa system was implemented in the simulation code, used in our research 
Group. The idea to apply the HNC method to produce highly optimized Jastrow terms 
of trial wave functions can be used in future Monte Carlo simulations of other physical 
systems. 

In the next part of the work, given in the Thesis, we presented a quantum Monte 
Carlo study of a bulk system of bosonic Rydberg atoms, that is alkali atoms with a single 
electron residing on a very high orbital, interacting through van der Waals interaction 
potential l/|r| 6 . A perspective utilization of Rydberg atomic clouds as quantum gates or for 
an observation of Bose-Einstein condensation, requires an utmost stable and predictable 
system, that often implies that the interaction is repulsive and independent on spacial 
directions, hence our simple model should be physically relevant. The asymptotic van 
der Waals C*6 / 1 1" | 6 is a common leading term of the interaction of Rydberg atoms, derived 
by virtue of the perturbation theory. One of the aims that we pursue in the study is 
to understand how much of the behavior of a real system comes from this leading term, 
and to which extent the l/|r| 6 interaction of the excited atoms in the cloud can describe 
crucial properties of experimentally relevant systems like mixtures of excited and ground 
state atoms. 

The zero-temperature properties of a system of Rydberg atoms, as presented in Sec- 
tion 5, are controlled by a unified dimensionless parameter prjj, which depends on the 
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density p and the interaction strength r$ = (mCe/h 2 ) 1 ^. The behaviors of distinct sys- 
tems with different parameters like particle mass and interaction strength constant, are 
therefore identical, if they possess the same value of prjj, for instance, the systems solidify 
and melt at the same dimensionless density. Making use of the diffusion Monte Carlo 
method, we found that for a system, governed by a model Hamiltonian 

i i<j 'i 3 I 

rewritten in a proper dimensionless form, the zero-temperature transition point happens 
in the density range 

p c rl = 3.9 ±0.2. 

The position of the phase transition is obtained by applying the double-tangent Maxwell 
construction technique implying that the pressure and the free energy stay unchanged 
along the melting curve. Possible types of the crystal packing, preferred by the system in 
its solid state, have been discussed. The calculations of the Madelung energy for different 
lattices suggest that face-centered cubic (fee) and hexagonal close-packed (hep) are the 
preferable packings with a slight advantage of fee, although the introduction of a quantum 
defect correction 8 in the interaction potential can give the advantage to hep. We stick 
to the fee lattice, which was used throughout our quantum Monte Carlo simulations as a 
preferable crystal formation in the solid phase. 

The finite temperature properties of the system were studied both with classical and 
quantum Monte Carlo methods. First of all we applied the classical Monte Carlo technique, 
based on the evolution of the system in accordance with its classical partition function. The 
estimates for the liquid-solid phase transition curve in case of high temperatures, where 
the classical approach is valid, have been found. The region of low temperatures, where the 
quantum description is required, was studied with the path integral Monte Carlo (PIMC) 
method. This technique allowed to find the region of the solid-gas phase transition as 
well as to localize the transition between normal fluid and superfluid phases. The results 
of the PIMC simulations completed the phase diagram of the system as a function of 
dimensionless density and dimensionless temperature, on the other hand confirming the 
correctness of the DMC and classical calculations, as the PIMC data came as a smooth 
transition between the two. 

We also present a discussion of a possible treatment of the crystallized Rydberg atom 
clouds as a model for perspective research of the supersolid phase. 

In Section 6 we performed extensive quantum Monte Carlo calculations of atomic para- 
hydrogen p-H.2 at zero temperature below its solidification curve in a range of densities. 
In this study our principal motivation was to understand better the properties of the 
metastable liquid/glass phase at low temperatures. In the limit of zero temperature we 
have used the DMC method, which is a very efficient tool to sample this metastable phase 
through the use of a trial wave function with a corresponding symmetry. The data provided 
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by diffusion Monte Carlo method suggests that the Bose-Einstein condensate fraction is 
subject to a large depletion compared to that of the stable liquid 4 He. 

In conclusion in this Thesis we successfully applied quantum Monte Carlo techniques for 
studying the quantum phase diagrams in a number of physical systems which are relevant 
to current and future experiments. Our research Group 

• learned to carry out the Ewald summation for different physical systems and to use 
it for solving the finite size correction problem; 

• studied and implemented in the code the modified periodic boundary conditions 
(truncated tetrahedron), which might be used to enhance the efficiency of our Monte 
Carlo calculations; 

• realized the usefullness the hypernetted chain (HNC) method for constructing opti- 
mized Jastrow terms of trial wave functions; 

• improved the understanding of the ways to study quantum phase transitions by virtue 
of quantum and classic Monte Carlo method. 

As a consequence of this work, and as a result of this Thesis, we published the following 
articles: 

1. O. N. Osychenko, G. E. Astrakharchik, Y. Lutsyshyn, Yu. E. Lozovik, and J. 
Boronat: "Phase diagram of Rydberg atoms with repulsive van der Waals interac- 
tion' 1 , Phys. Rev. A 84, 063621 (2011). 

2. O. N. Osychenko, G. E. Astrakharchik, and J. Boronat: "Ewald method for polytropic 
potentials in arbitrary dimensionality" , Mol. Phys. 110, 4, 227-247 (2012). 

3. O. N. Osychenko, G. E. Astrakharchik, F. Mazzanti, and J. Boronat: "Zero-temperature 
phase diagram of Yukawa bosons", Phys. Rev. A 85, 063604 (2012). 



4. 



O. N. Osychenko, R. Rota, and J. Boronat: 11 Superfluidity of metastable glassy bulk 
para-hydrogen at low temperature" , Phys. Rev. B 85, 224513 (2012). 



Appendix A 

Ewald method for polytropic 

potentials 



We prove that the sums 5_ + and S ++ (3.46-3.47) vanish on average, allowing to calculate 
the potential energy over the negatively charged particles' positions only. 

• First, let us show that the integral of ip over the cell is zero. Since the distances are 
in the units of L, consider the cubic cell fl = (x,y, z) G [—1/2, 1/2] 3 , that yields 

/ ij(r)dr = J 1 + J 2 + C 1 (A.l) 
Jn 

Ji = I dr^i?(n,r) (A.2) 

J 2 = J n dr J2 K ( m > r ) • ( A - 3 ) 



where 



It can be easily seen, that the second integral J 2 is zero 



J 2 = x(m, r) / cos(27rmr) dr 



r^o (2itym x m y m z 

As far as the integral J\ is concerned, we can notice that the regions Q'(n) = r + n, 
where r G Q, n G 1? are the same cubic unit cells, displaced by an integer vector, 
thus covering all the coordinate space with only zero-measure intersections. It means 
that the summation of the integrals in (A.2) over the cell Q can be substituted by 
the integration over the whole coordinate space, 

k „,2 1 



r x — ^ f r(f, « k + n\ 2 ) 

J x =Y R(n, r)dr = Y / V2 ' ,' 1 



r(f)p fc (fc-3)r 



and thus the whole integral (A.l) is equal to zero. 
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• Consider two species of the particles: negative charges g» on positions and a 
positively charged and uniformly distributed background with a total charge q+N + = 
—qiNi, ensuring the neutrality of the cell. Let us demonstrate that SL + is equal to 
zero, when the number of background charges tends to infinity. In this case the sum 
(3.46) for S |_ may be rewritten as an integral over the background charges' positions 

S-+ = YjQi / ^( r p ~ r i) adr P = / ^(r-)crdr , (A.6) 

where we did the change of variables r = r p — r.;. The regions Q and Qi refer 
to the original simulation cell and the same cell, moved by the vector 7*3, and a 
stands for the background charge density o = —qiNi/V(Q). It is clear that every 
vector r = (2, y, z) 6 f2,- can be displaced into the cell Q by the corresponding shift 
r = (x, y, z) — (x — aL, y — bL, z — cL) G Q with integers a, b, c. The Jacobian 
J of the change of variables r — > r is obviously 1. On the other hand, due to the 
periodicity of if), 

^(r)=^(f), (A.7) 

and f runs over the whole region Q due to the conservation of the volume with J = 1. 
Finally, Eq. (A.6) can be written as 

/ tp(r)adr = 0. (A.8) 

In the similar manner, the interaction between the charges of the background S ++ 
in the limit — > 00 is given by the double integral 

S++ = l [ dn I dr 2 i>( ri - r 2 )<j 2 = , (A.9) 
2 Jn Jn 

since J n ij){r\ — r 2 ) dr 2 = 0, following the same arguments as for the case of S- + . 



Appendix B 

Truncated octahedron boundary 

conditions 



Apart from the standard cubic periodic boundary conditions we expored a possibility 
to apply the truncated octahedron periodic boundary conditions. This kind of periodic 
conditions has an advantage to reduce effects of anisotropy as well as enhance the overall 
efficiency of the simulation. 

The truncated octahedron simulation cell, represented by a cube with removed angles, 
possesses a periodicity with the steps {L} p = (±L/2, ±L/2, ±L/2), which means that 
only a crystal formation of a cubic class (simple cubic, bcc, fee) with a composition, 
corresponding to even values of L/a (a = (N ce u/ p) 1 ^ 3 stands for a size of the elementary 
cell) is commensurate with the cell. It means that the simulation in periodic boundary 
conditions is only possible for a certain reduced set of particles in the box. In practice 
it is equivalent to the exclusion of the crystal configurations with odd values of L/a; for 
instance, in case of fee lattice one has to omit N p = 2(2i + l) 3 = 2, 64, 250, 686.... 

The situation regarding hep formation is similar, with the only difference that as it 
always happens in the case of hep the simulation should be performed in a cell with unequal 
size lengths L x , L y , L z . If the truncation of the interaction (minimum image convention) 
is to be applied, one takes the distance from the center to the nearest plane, that is the 
plane with the equation X/L x + Y '/ L y + Z / L z = 3/4. The spherical cut-off radius in this 
case is equal to 

3/4 



cut 



l/Ll + l/Ll + l/Ll 



(b.i; 



The values of L x , L y , L z are chosen such that the cell is commensurate with the lattice. 
We also require that an arbitrary periodic translation r — > r + {L} brings a lattice site 
again to a site. This condition results in a demand of having all numbers L x /a x , L y /a y , 
L z /a z even, that again reduces the variety of options for a number of particles in a cell to 
N p = 32,256,864,.... 

The number of particles N p in a simulation in a truncated octahedron cell corresponds 
to a cut-off radius 

r c =^(2N p /p)^ (B.2) 

which is a distance to the nearest plane of the cell. The same distance in the cubic p.b.c. is 
obtained by 3\/3/4A^ ps 1.3 or, as diffusion Monte Carlo technique is 0(N p ), by about 68% 
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of additional calculation time. This effect is hindered by a more complex implementation 
of a particle motion, but in practice the overall efficiency gap is never below 30%. 

As it was mentioned above, a truncated octahedron cell replicated periodically with 
the displacements from the set {L} fills entirely the coordinate space. This makes the 
application of the Ewald technique theoretically possible, once the original potential energy 
sum is rewritten in a suitable form, compatible with the standard x, y, z-axis periodicity, 
required by the Ewald method. Indeed, as one can easily notice, the total set of cell images 
can be classified into two groups: the images, produced by an even number steps with a 
resulting cell displacement (iL, jL, kL) (i, j, k are arbitrary integers) and by an odd number 
of step with a displacement (iL, jL, kL) + (L/2, L/2, L/2) . If one thinks of a simulation cell 
as a compound of the original cell fl and the displaced cell fid = {L/2, L/2, L/2) + fl, the 
total space will be covered again by its replications along the axis x, y, z with the period 
L. This allows to apply the Ewald summation technique directly to this complex cell with 
the number of particles equal to 2N p , but due to the similarity of fl and f^, the sums may 
be simplified in the following way: 



2N„ 



5 



-Y. 

2 ^ 

* n 



E 

n 

E 

TO 

E 



E + E +E E + E E 

V=Wj i,j=N p +l;i^j i=lj=N+l i=N p +lj=l 



E +E E 
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i.i < 



£ Rim) + £ R( rij + L/2) 



£ (R(r l3 ) + R(r {j + L/2)) + J2 R(L/2) 



R{ rij ) 



(B.3) 



where we place the points 1,...,N P to the cell Q and their corresponding images 1 + 
Np, 2Np to fid- I n a similar line of thought we can represent the Fourier-transformed 
part of the Ewald sum: 



^ 2N P 



E E Kfa) = E ( E Wry) + K( nj + L/2))) + K(L/2)) . 



(B.4) 



8=1 



The constants £, C\, C 2 (refer to Eqs. (3.36,3.38,3.39)) stay clearly unchanged, as they 
characterize the self-image interactions, not biased by a specific cell geometry. Notice that 
the number of particles, appearing as a factor of £ in the expression for the total energy is 



now 2N P . 
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